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Preface 

This  thesis  reports  two  studies  concerned  with  the 
feasibility  of  plasma  containment  In  toroidal  magnetic 
fields.  Lt.  Col.  R.  C.  Wingerson  of  the  Plasma  Physics 
Research  Laboratory  of  the  Aerospace  Research  Laboratories 
sponsored  this  research.  Dr.  Wingerson’s  Interest  In  mag¬ 
netic  containment  schemes  Is  reflected  by  the  current 
experimental  and  theoretical  work  of  the  Plasma  Physics 
Research  Laboratory,  to  which  this  report  might  hopefully 
have  some  application.  Much  of  the  theoretical  work  upon 
which  this  thesis  Is  based  Is  Dr.  Wingerson’s  own. 

The  first  part  of  the  study  Is  an  Investigation  of  the 
magnetic  field  of  an  Infinite  system  of  thin  loops  of  cur¬ 
rent.  The  second  part  develops  the  equations  describing 
an  axlsymmetrlc  system  of  charged  particles  In  crossed 
electric  and  magnetic  fields. 

I  would  like  to  thank  Dr.  Wingerson  lor  his  very 
generously  given  time  and  effort.  I  am  also  Indebted  to 
the  Applied  Mathematics  Research  Laboratory  and  the  Digital 
Computation  Division  of  the  Aeronautical  Systems  Division  at 
Wright- Patters on  Air  Force  Base.  My  wife  has  contributed  a 
great  deal  of  encouragement  and  understanding. 

Richard  D.  Franklin 
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Abstract 

A  toroidal  magnetic  plasma  containment  configuration 
is  proposed  wherein  coils  are  centered  on  a  closed  helix. 

It  is  possible  to  optimize  the  qualities  of  the  field  by  ad¬ 
justing  the  coll  configuration.  The  equations  describing 
an  equilibrium  plasma  in  an  axlsymmetric  system  are  derived 
in  terms  of  the  particle  density  distributions,  radial  and 
azimuthal  drift  velocities,  and  the  electric  and  magnetic 
field  strengths  necessary  to  maintain  equilibrium. 
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CONTAINING  CONTROLLED  FUSION  REACTIONS 
WITH  CROSSED  ELECTRIC  AND  .1AGNETIC  FIELDS 

I .  Introduction 

Dating  back  to  1944  the  possibility  of  controlled 
fusion  has  held  the  Interest  of  a  large  segment  of  the 
world  scientific  community .  The  basic  principles  of  con¬ 
trolled  thermonuclear  reactions.  If  not  already  known,  ware 
laid  down  by  scientists  at  the  Los  Alamos  Scientific  Labo¬ 
ratory  as  early  as  1946.  Among  these  were  Fermi,  Teller, 
Tuck,  and  Wilson.  Since  that  time,  however,  no  controlled 
fusion  reactions  of  extended  duration  or  energy  have  been 
possible. 

Confinement 

One  of  the  major  obstacles  to  controlled  fusion  Is 
still  the  need  for  a  method  of  plasma  confinement.  Almost 
all  currently  proposed  methods  can  be  divided  Into  three 
types  of  devices i  the  pinch  (Ref  2t22),  the  magnetic 
mirror  (Ref  2«51),  and  the  torus  (Ref  2i33).  All  of  these 
utilize  magnetic  fields  In  attempting  to  trap  the  charged 
particles  composing  a  plasma  Into  helical  orbits  about 
field  lines.  The  objective  of  all  methods  Is  to  contain 
the  colliding  particles  long  enough  for  a  copious  number  of 
nuclear  reactions  to  occur. 
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A  high-temperature  plasma  Is  much  like  an  ordinary  gas 
in  that  it  tends  to  diffuse  because  of  interparticle  colli¬ 
sions.  Magnetic  containment  devices  are  designed  to  exert 
containing  pressures  on  the  expanding  plasma  and  limit 
diffusion.  It  can  be  shown  (Ref  5«275)  that  the  rate  of 
plasma  diffusion  across  a  straight,  uniform  magnetic  field 
is  inversely  proportional  to  the  square  of  the  magnetic 
field  strength.  This  is  important  because  It  shows  that  If 
stable  confinement  is  possible  diffusion  can  be  reduced  by 
raising  field  strength.  Thus,  If  the  density  of  the  plasma 
can  be  raised  high  enough,  thermonuclear  reactions  can 
occur  at  a  self-sustaining  rate. 

The  toroidal  magnetic  field  is  attractive  for  plasma 
confinement  because  it  presumably  does  away  with  the  end 
effects  of  the  pinch  and  the  magnetic  mirror.  The  toroid 
appears  to  the  plasma  as  an  endless  tube.  In  this  type  of 
system  the  field  is  generated  with  current  windings  around 
a  toroidal  shape.  The  field  lines  are  roughly  circles 
inside  the  torus.  It  is  easily  seen,  however,  that  the 
field  within  the  torus  is  stronger  at  the  inside  circum¬ 
ference  than  at  the  outside  of  the  torus  because  the 
windings  are  closer  together. 

Several  toroidal  devices  have  been  proposed  and  built 
to  correot  this  non-uniformity.  Among  them  are  the  Scyllac 
(Ref  13«5^3)*  the  Wisconsin  (Ref  llilll5),  and  the  most 
famous,  the  Stellarator  (Ref  2i37)« 
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In  1951  Spltzer  proposed  a  device  oalled  the 
"Stellarator"  which  Introduced  a  rotation  of  field  lines 
within  a  twisted  torus.  This  phenomenon  cancelled  charged 
particle  drifts  transverse  to  the  field.  Later  he  plaoed 
helical  windings  around  an  ordinary  torus  which  accomplished 
the  same  thing. 

Part  of  this  thesis  Is  a  study  of  the  field  of  a 
toroidal  device  proposed  by  Lt.  Col.  R.  C.  Wlngerson  of  the 
Plasma  Physics  Research  Laboratory  of  the  Aerospaoe  Research 
Laboratories.  This  device  Is  a  series  of  solenoid  colls 
with  centers  roughly  conoentrlc  on  a  circle.  The  oenters  of 
the  colls  are  displaced  off  the  major  circumference  of  a 
torus  by  a  small  distance  In  such  a  way  as  to  leave  an  equal 
angle  between  the  directions  of  displacement  of  successive 
colls.  Pig.  1  Is  a  schematic  representation  of  this  system. 
The  parameters  of  the  c onf  1  gurat  1  on ,  such  as  the  number  of 
colls,  the  spacing  between  colls,  the  displacements  off 
axis,  and  the  angle  between  displacement,  are  the  variables 
used  to  manipulate  the  field. 

If  only  a  short  segment  of  the  arc  of  this  toroid  Is 
considered,  the  system  can  be  approximated  by  an  Infinite ly- 
long  straight  system  of  colls.  The  colls  are  then  plane- 
parallel. 

Chapter  tl  of  this  paper  describes  the  methods  and 
procedures  for  the  study  of  the  field  of  this  Infinite 
straight  system,  and  lt  reports  the  findings.  The  analysis 
of  the  field  Is  directed  mainly  towards  determining  the 
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Fig.  1.  Schematic  of  Toroid  of  Displaced  Colls 
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transverse  uniformity  of  the  field  strength,  the  linearity 
of  the  field  lines  and  the  stellarator-type  rotations  of  the 
field  lines. 

The  purpose  of  this  study  is  to  determine  the  feasibil¬ 
ity  of  designing  a  theoretically  useful  containment  device 
by  adjusting  the  parameters  of  the  coil  configuration.  If 
the  characteristics  of  the  field  can  oe  controlled  by  vary¬ 
ing  the  spatial  parameters  of  the  coil  configuration,  we 
may  have  a  feasible  way  of  designing  magnetlo  fields  with 
very  specific  properties. 

Plasma  in  Equilibrium 

The  second  part  of  this  study  develops  the  hydro- 
dynamic  equations  of  a  system  of  two  types  of  particles, 
electrons  and  positive  ions,  in  equilibrium  with  an  axially 
directed  magnetic  field.  Only  the  time -Independent  oase  is 
dealt  with.  First  a  generalized  distribution  function  is 
derived  in  which  particle  energy  and  canonical  angular 
momentum  are  conserved.  The  resulting  function  is  essen¬ 
tially  Maxwellian,  but  an  azlmuthally  symmetric  rotation  of 
the  single  distribution  is  predicted.  It  is  necessary  to 
derive  an  average  volumetric  foroe  whioh  arises  from  colli¬ 
sions  between  unlike  partloles  in  order  to  aooount  for  the 
interaction  of  the  superimposed  distributions. 

Tho  steady-state  Boltzmann  equation  is  relied  upon  to 
provide  the  hydrodynamic  equations  for  each  distribution. 
Exoept  for  the  inclusion  of  a  foroe  term  describing  changes 
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In  the  distributions  due  to  collisions,  the  distributions 
are  treated  entirely  separately.  One  of  the  objectives  of 
the  study  is  to  test  theoretically  the  effects  on  the 
equilibrium  situation  of  great  differences  between  the 
electron  and  ion  kinetic  temperatures. 

Another  prime  objective  of  the  study  is  to  investigate 
the  feasibility  of  prescribing  radial  particle  number 
density  profiles  by  the  injection  of  charged  particles.  It 
is  assumed  that  electrons  and  ions  can  be  separately  depos¬ 
ited  in  the  plasma  at  any  radial  position  and  at  any  time 
rate  desired. 

If  moments  of  the  Boltzmann  equation  are  taken  in 
cylindrical  coordinates,  all  the  variables  describing  each 
distribution  become  averages  which  are  functions  of  the 
only  Independent  variable,  radius.  The  two  non-zero 
average  velocities  are  the  azimuthal  and  radial  velocities. 
Density  is  a  third  variable  describing  each  distribution. 
The  moments  of  the  Boltzmann  equation  yield  one  continuity 
equation  and  two  momentum  transfer  equations  for  each 
distribution.  These  six  equations  are  all  first  order 
differential  equations. 

Two  of  Maxwell’s  equations  apply  in  the  steady-state. 
Current  densities  and  charge  density  are  finite  but  are 
taken  as  looal  averages.  The  magnetic  field  is  assumed 
axially  directed,  and  its  magnitude  on  axis  is  to  be  se- 
leoted.  The  eleotrio  field  is  assumed  entirely  radial  in 
direction,  and  no  external  electrlo  field  is  utilized.  In 
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the  cylindrical  geometry  then  Maxwell's  equations  supply 
two  first  order  differential  equations. 

Thus  a  total  of  eight  simultaneous  differential  equa¬ 
tions  are  available.  They  are  in  terms  of  the  two  densities, 
the  pairs  of  radial  and  azimuthal  average  velocities,  and 
the  electric  and  magnetic  fields.  Thus  there  are  in  prin¬ 
ciple  enough  differential  equations  to  solve  simultaneously 
for  the  radial  profiles  of  the  distribution  variables  and 
the  electric  and  magnetic  fields.  Since  the  distribution 
function  is  based  on  an  equilibrium  system,  these  profiles 
should  describe  an  equilibrium  plasma. 

Procedures  are  to  be  developed  for  studying  the  effects 
on  such  an  equilibrium  of  the  injection  of  particles,  and  of 
the  imposition  of  temperatures,  and  the  axial  magnetic  field. 
Involved  is  the  development  of  numerical  integration  methods 
for  the  simultaneous  solution  of  the  equations.  A  careful 
analysis  is  made  of  those  characteristics  of  the  hydro- 
dynamic  equations  which  determine  the  applicability  of 
numerical  methods  to  their  solution. 

The  type  of  plasma  of  interest  in  this  study  is  the 

high-temperature,  high-density  thermonuclear  ionized  gas. 

It  Is  current  opinion  that  temperatures  on  the  order  of 

10^  -  10^°K  and  particle  densities  on  the  order  of  lO1^ 

-3 

cm  are  required  to  sustain  fusion  reactions  in  the  plasma 
environment.  In  Chapter  IV  some  numerical  examples  are 
carried  out  to  examine  some  of  the  predicted  velocities  and 
electric  fields  which  correspond  to  this  type  of  plasma. 
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Much  of  the  theoretical  work  leading  to  the  derivation 
of  these  equations  was  done  by  Dr.  Wlngerson  In  unpublished 
papers.  He  and  the  scientists  of  the  Plasma  Physics 
Research  Laboratory  are  at  this  time  actively  engaged  In 
the  design  and  testing  of  axlsymmetrlc  plasma  systems.  It 
Is  hoped  that  the  efforts  reported  in  this  paper  can  be 
applied  to  the  design  of  experimental  plasma  containment 
devices . 
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II.  Magnetic  Fie Id  of  the  Dlsplaced-Coll  C onf lgurat 1 on 

Introduction 

This  chapter  Is  a  report  of  a  study  on  the  applica¬ 
bility  of  a  proposed  plasma  containment  device,  dubbed  the 
"dlsplaced-coll  configuration."  The  major  objective  of  the 
study  was  to  determine  how  the  characteristics  of  the  field 
of  a  device  similar  to  the  system  of  colls  shown  In  Fig.  1 
depend  on  the  geometric  arrangement  of  the  separate  oolls. 
Some  of  these  Important  characteristics  are  mentioned  In 
Chapter  Ij  of  Interest  are  the  gradients  of  field  strength 
through  the  field,  the  direction  of  the  field,  and  the 
stellarator  rotations  of  field  lines  In  the  field.  It  was 
not  the  objective  of  the  study  to  determine  what  field 
characteristics  are  necessary  or  desirable  for  plasma  con¬ 
tainment.  Nor  was  the  objective  to  design  a  magnetic  field 
using  the  dlsplaced-coll  configuration.  This  study  was 
made  to  ascertain  the  feasibility  of  controlling  the  major 
characteristics  of  a  toroidal  magnetic  field  by  varying  the 
geometric  parameters  of  the  dlsplaced-coll  configuration. 

First,  the  geometry  of  a  straight  coll  configuration 
which  approximates  the  configuration  of  Fig.  1  Is  defined. 
Notation  of  variables  Is  established.  An  orthogonal  coor¬ 
dinate  system  Is  also  defined  relative  to  the  coll  system. 

The  equations  for  the  vector  components  of  B  are 
derived,  and  computer  programs  are  described  which  calcu¬ 
late  the  magnitude  and  direction  of  B.  Procedures  are 
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developed  to  systematically  Investigate  field  variables  as 
functions  of  field  position  and  the  geometric  parameters. 

The  method  of  generating  mathematically  the  field  lines 
of  the  magnetic  field  Is  explained.  Two  characteristics  of 
Interest  In  the  pattern  of  field  lines  are  Identified,  and 
the  procedures  for  studying  each  are  explained. 

The  section  covering  results  reports  the  variations  of 
the  quantities  of  Interest  through  the  field,  and  their 
functional  relationships  to  the  geometric  parameters.  The 
feasibility  of  designing  fields  with  very  particular 
characteristics  Is  ascertained. 

Coll  C onf lgurat 1  on 

Let  us  consider  as  an  approximation  to  the  torus  of 
real  multi pie -turn  colls  an  infinite  straight  system  of 
thin,  single-loop  colls.  The  approximation  Is  valid  If 
only  a  small  segment  of  arc  of  the  torus  of  solenoids  Is 
considered. 

The  magnetic  field  of  a  single  loop  of  current  is  very 
similar  to  that  of  a  finite  solenoid.  The  strength  of  the 
B  field  varies  similarly  In  either  as  a  function  of  posi¬ 
tion,  and  the  shapes  of  the  fields  of  the  two  are 
characteristically  the  same. 

The  equations  of  the  transverse  and  social  components 
of  B  for  an  Inf lnltely- thin  loop  of  01  r.-ent  are  presented 
In  Appendix  A.  Techniques  are  der?  .  a  there  for  computing 
numerloal  values  for  these  components. 
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Fig,  2  Is  a  schematic  representation  of  the  scheme 
for  constructing  the  infinite  straight  system  of  colls.  The 
axial  Interval  between  colls  Is  designated  "space."  The 
displacement  of  the  coll  centers  off  the  Z  axis  Is  "dlsp." 
These  equal  displacements  are  made  such  that  there  Is  a 
constant  angle  "alpha"  between  the  directions  of  dlsp  of 
neighboring  colls.  The  arbitrary  convention  is  made  that 
the  slgr.  of  alpha  Is  positive  as  measured  from  one  coll 
center  to  the  next  when  proceeding  from  coil  to  ooll  In  the 
positive  axial  direction.  The  radius  of  all  colls  Is  "a." 

The  result  of  these  manipulations  Is  a  configuration 
of  plane-parallel  loops  whose  centers  lie  on  a  uniform  he¬ 
lix.  The  pitch  length  of  the  helix  Is  the  product  of  space 
and  alpha,  and  its  radius  is  dlsp.  The  axis  of  the  helix  Is 
chosen  as  the  Z  axis  of  am  orthogonal  coordinate  system. 

The  X  axis  of  the  system  Is  arbitrarily  located  through  the 
center  of  one  of  the  colls. 

Let  the  radii  of  the  colls  be  normalized  to  unity,  and 
let  all  distances  henceforth  be  measured  in  units  of  coll 
radius.  Choose  the  currents  In  the  colls  to  be  equal,  In 
the  same  direction,  and  of  magnitude  such  that 


*pl  -  \ 

2_tt 


(1) 


This  concludes  the  definition  of  the  dlsplaced-coil 
configuration.  Its  geometric  parameters  are  space,  dlsp, 
and  alpha.  This  chapter  Is  primarily  concerned  with  the 
dependence  of  the  magnetic  field  on  these  three  variables. 
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Objectives  of  Study 

This  section  describes  more  fully  the  areas  to  be 
Investigated  In  this  study.  It  also  attempts  to  supplement 
Chapter  I  In  Its  Justification  of  the  study. 

Characteristics  of  B.  Again,  the  properties  of  B  most 
critical  to  containment  feasibility  are  uniformity  of 
strength  and  the  ratio  of  transverse  to  axial  component 
strengths.  This  study  attempts  first  to  discover  In  general 
the  orders  of  magnitude  of  these  properties.  Then  relations 
between  ;heo  and  the  three  conf lgurat ion  parameters  are 
sought.  The  variations  of  the  two  properties,  as  functions 
of  both  the  radius  out  from  the  axis  and  position  along  the 
axis,  are  also  Investigated. 

Field  Line  Path3 .  The  magnetic  field  lines  surround¬ 
ing  a  3lngle  thin  loop  are  closed  loop3  In  plane3 
perpendicular  to  the  plane  of  the  loop.  Field  strengths 
vary  along  the  paths  of  these  lines,  but  the  magnitudes  are 
gymnetrlc  about  the  plane  of  the  loop. 

When  many  single  loops  are  arranged  with  their  centers 
on  a  helix  the  field  lines  of  the  field  can  .no  longer  be 
planar.  The  field  Is  expected  to  exhibit  a  helical  twist 
In  the  region  Inside  all  colls. 

Questions  Immediately  come  to  mind.  •  Do  the  field 
lines  encircle  the  axis  of  the  system  as  does  the  helix 
describe!  by  the  coll  centers,  or  do  they  remain  In  the 
same  quadrant  of  the  system?  Do  subharmo.olcs  of  this  field 
line  twist  exist  such  that  the  entire  field  line  pattern 
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rotates  as  a  solid  body?  No  matter  how  either  of  these 
rotations  appear  their  frequencies  and  the  radii  of  the 
resulting  helices  should  be  determined. 

Methods  and  Procedures 

This  section  outlines  the  computational  methods  and 
procedures  employed  In  studying  the  magnetic  field. 
Appendices  are  used  to  support  this  explanation.  Appendix  A 
Is  the  derivation  of  the  numerical  approximations  for  the 
transverse  and  axial  vector  components  of  B.  Appendix  B 
describes  a  computer  subroutine  called  "Field''  which  com¬ 
putes  the  direction  cosines  and  strength  of  B  at  any  given 
point  in  the  field.  The  program  "Tracer”  Is  explained  In 
Appendix  Ci  It  Is  designed  to  Integrate  the  paths  of  the 
field  lines. 

Ranges  of  Variables.  The  volume  of  Interest  In  the 
magnetic  field,  and  the  ranges  over  which  the  geometric 
parameters  are  to  be  studied  need  specification.  Because 
of  the  undesirable  behavior  of  the  field  variables  near  the 
conductors,  the  volume  over  which  the  field  Is  expected  to 
be  at  all  well-behaved  Is  limited  to  be  entirely  within  all 
conductors.  The  size  of  dlsp  obviously  determines  the  size 
of  the  volume  satisfying  this  condition. 

The  ranges  of  the  geometric  parameters  to  be  considered 
in  this  study  are  as  follows i 

spacel  0.1  —  1.0 

dispt  0.  —  0.2 

alpha i  0  —  lBo  degrees 

14 


GS P/PH/6 9- 7 


Limits  are  Imposed  on  space  and  dlsp  only  for  definiteness, 
and  are  not  meant  to  be  exclusive.  They  can  be  expanded  if 
they  do  not  seem  to  Include  all  the  interesting  cases.  By 
symmetry  the  range  of  alpha  covers  also  the  cases  for  alpha 
greater  than  180  degrees.  The  above  intervals  are  sub¬ 
divided  sufficiently  to  reveal  the  important  functional 
relationships . 

Bm  and  RTA.  Let  "Bm"  refer  to  the  magnitude  of  B. 
Appendix  A  presents  the  equations  for  the  axial  (Bz)  compo¬ 
nent  and  the  total  component  transverse  to  the  axis  (Br). 

Let  "RTA”  stand  for  the  ratio  of  Br  and  Bz.  Field  is  a  sub¬ 
routine  which  can  deliver  Bm  and  RTA  at  any  location. 

The  procedure  is  to  compute  both  these  terms  for  any 
given  set  of  space,  dlsp,  and  alpha  at  several  different 
radii  within  all  loops  and  also  at  points  along  the  axis 
between  any  two  loops.  The  computations  for  varying  radius 
are  conducted  halfway  between  two  colls  to  minimize  anom¬ 
alies  due  to  nearby  conductors. 

The  analysis  of  results  consists  first  of  finding  how 
RTA  varies  with  space,  dlsp,  alpha,  radius  from  the  axis, 
and  position  along  the  axis.  Then  the  rates  of  change  of 
Bm  with  radius  and  axial  position  are  determined  as  gener¬ 
alized  functions.  Finally,  these  rates  of  change ,  <)Bm/<)r 
and  <)Bm/iz,  are  investigated  as  functions  of  space,  dlsp, 
and  alpha.  Since  the  normalized  Bm  is  a  strong  function  of 
space,  it  is  necessary  to  calculate  dBm/dr  and  dBm/dz  as 
relative  rates  of  change. 
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Field  Line  Pattern.  We  have  mentioned  that  the  follow¬ 
ing  two  phenomena  are  probably  of  Interest i  the  helical 
twist  of  the  lines  caused  by  the  helical  coll  conf lgurat 1  on , 
and  the  possible  slow  rotations  of  the  entire  field  line 
pattern  about  the  axis.  Let  us  call  the  helical  twist  a 
first  order  occurrence.  This  twist  is  studied  as  a  function 
of  radial  position  in  the  field  and  of  the  coll  geometry. 

Field  lines  are  traced  over  one  repetition  or  cycle  of 
coll  displacement.  The  tracings  are  begun  from  points  at 
various  radii  In  the  plane  of  an  arbitrary  coll. 

A  convenient  measure  of  the  resulting  helices  Is  the 
helical  radius.  It  Is  expected  that  the  pitch  length  of 
the  first  order  twist  is  the  same  as  that  of  the  coll 
pattern.  That  Is,  the  pitch  length  of  the  lines  should  be 
the  same  as  that  of  the  coll  center  helix.  This  must  be 
verified,  however. 

It  Is  necessary  to  determine  If  the  second  type  of 
rotation  Is  Important  to  the  shape  of  the  field.  These 
rotational  drifts  become  evident  only  over  several  cycles 
of  the  first  order  twist.  Therefore,  a  few  lines  beginning 
at  different  radii  are  traced  over  one  hundred  first  order 
cycles  for  several  sets  of  space,  dlsp,  and  alpha.  If  the 
frequency  of  the  rotation  of  any  line  around  the  axis  Is 
such  that  the  axis  is  encircled  In  less  than  one  hundred 
oyoles,  the  second  order  rotation  may  be  important.  A  more 
detailed  analysis  of  the  frequencies  of  these  rotations  Is 
then  to  be  undertaken. 
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The  field  lines  are  generated  by  successively  Integrat¬ 
ing,  point  ty  point,  the  direction  cosines  of  B  over 
Incremental  distances.'  Appendix  C  describes  and  gives 
operating  Instructions  for  a  program  which  Integrates  lines, 
beginning  at  selected  X-Y  coordinates,  over  selected  dis¬ 
tances. 

The  program  supplies  a  record  of  the  three  coordinates 
of  points  along  the  field  line.  It  Is  also  written  to  be 
amenable  to  plotting  subroutines  so  that  the  lines  can  be 
traced  In  machine  drawings. 

The  diameters  of  the  first  order  helices  are  easily 
calculable  from  the  printed  records  of  the  X  and  Y  coordi¬ 
nates  of  a  line.  We  need  merely  subtract,  say,  the  minimum 
X  coordinate  from  the  maximum.  Tracer  supplies  these  coor¬ 
dinates  ;o*four  places. 

The  pitch  length  of  the  lines  Is  the  same  as  that  of 
the  coil-center  helix  If  the  field  lines  cross  the  planes 
of  concentric  colls  at  the  same  X  and  Y  coordinates. 

Appendix  C  details  the  maximum  errors  lr.  these  coordinates 
which  ca.i  arise  cut  of  the  numerical  Integration  methods 
utilized 

To  measure  the  frequencies  of  the  drifts,  the  numerical 
listing  of  point  coordinates  Is  studied  to  find  the  differ¬ 
ences  in  both  the  X  and  Y  coordinates  of  a  line  as  it  passes 
through  the  planes  of  colls  one  hundred  first  order  cycles 
apart.  The  resulting  differences  give  a  good  average 
measure  cf  the  magnitude  ar.i  direction  of  the  drift  of  a 
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part 1 cular  line. 

Results 

This  section  Is  a  report  of  how  the  characteristics  of 
the  field  of  the  dlsplaced-coll  configuration  vary  with 
space,  dlsp,  and  alpha.  The  dependence  of  RTA  and  the 
uniformity  of  Bm  on  these  parameters  Is  presented  graphical¬ 
ly  in  Figs.  3  and  4.  These  curves  quantitatively  depict 
RTA  and  the  relative  gradient  with  radius  of  Bm  on  axis, 
halfway  between  any  two  colls.  However,  curves  were  drawn 
for  other  locations  In  the  field,  and  their  shapes  are 
very  similar  to  the  curves  In  Figs.  3  and  4.  The  general 
dependence  of  RTA  and  Bm  on  position  in  the  field  Is  qual¬ 
itatively  described  below.  Also,  sample  computations  of 
RTA  and  dB/dr  for  values  of  space,  dlsp,  and  alpha  not 
Included  In  the  ranges  previously  set  down  indicated  no 
unexpected  behavior. 

RTA.  Fig.  3  contains  graphs  of  RTA  on  axis  as  a 
function  of  each  of  the  configuration  parameters  alone, 
holding  the  other  two  parameters  fixed.  RTA  Is  a  strong 
function  of  all  three  geometric  parameters  and  of  radial 
position  In  the  field.  The  almost  linear  function  of  RTA 
with  dlsp  Is  predictable.  Increasing  dlsp  accentuates  the 
helical  twist  of  the  field  lines,  and  hence  the  transverse 
components  of  B,  of  which  RTA  is  a  measure.  RTA  increases 
exponentially  with  radius  from  the  axis.  The  rate  of  In¬ 
crease  with  radius  depends  on  the  specific  combination  of 
spaoe,  dlsp,  and  alpha. 
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and  Alpha 
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Fig.  4,  dB/dr  versus  Space,  Dlsp,  and  Alpha 
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Bm.  The  mean  magnitude  B  is  inversely  proportional  to 
space,  and  is  not  strongly  dependent  on  disp  and  alpha. 

The  dependence  on  space  is  expeoted  slnoe  in  general  the 
magnetic  field  of  a  coil  varies  Inversely  with  distance 
from  the  center  of  the  coil. 

Uniformity  of  Bm.  Consider  first  the  relative  quantity 
dB/dr  on  axis  as  a  function  of  space,  disp,  and  alpha  as 
depicted  in  Pig.  4.  Note  that  the  curves  of  dB/dr  are  very 
much  like  the  RTA  curves.  Bm  does  not  vary  azlmuthally 
near  the  axis,  that  is,  well  within  all  the  colls.  With 
increasing  radius,  as  the  coils  themselves  are  approached, 
the  field  is  dominated  by  the  nearest  colls.  It  is  diffi¬ 
cult  to  describe  in  general  the  field  near  the  coils.  Bm 
is  quite  uniform  along  the  axis.  As  a  quantitative  exam¬ 
ple,  Bm  varies  along  the  axis  by  less  than  three  percent 
for  space  equal  to  unity. 

Field  Lines .  Fig.  5  shows  the  paths  of  ten  field 
lines  traced  over  two  cycles  of  coll  displacement.  Fig.  6 
is  a  listing  of  coordinates  along  the  line  nearest  the  axis 
for  two  cycles.  The  coll  configuration  for  these  lines  is 
defined  as  follows i 

s  pace  «  .  6 
disp  -  .2 
alpha  >■  60° 

The  pitch  length  of  the  line  twist  is  the  same  as  that 
of  the  coil  oenter  helix.  The  radius  of  the  helical  line 
on  axis  is  .07.  The  radii  of  the  helices  described  by 
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SPACE  BETWEEN  CfllCS  3.0. 
DISPLACEMENT  Of  COIL  CENTER  0.20 
ANGLE  BETWEEN  CENTERS  fci.J  DEGREES 
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Pig.  6.  Listing  of  Field  Line  Coordinates 
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lines  farther  off  the  axis  increase  exponentially  from  this 
value.  This  behavior  explains  the  dependence  of  RTA  on 
radius  from  the  axis.  The  pattern  of  the  first  order  twist 
shown  In  Fig.  5  Is  symmetric  about  the  axis. 

From  Fig.  5  It  appears  that  the  field  lines  do  not 
tend  to  rotate  about  the  Z  axis  of  the  system.  Actually 
there  Is  a  slight  drift  of  these  lines  about  the  axis. 
Selected  lines  were  integrated  over  one  hundred  first  order 
cycles  for  various  combinations  of  space,  dlsp,  and  alpha. 
The  transverse  coordinates  of  the  lines  did  not  return  to 
their  starting  values  after  each  first  order  revolution. 

The  indicated  drifts  over  one  hundred  first  order  cycles 
were  extrapolated  to  compute  how  many  first  order  cycles 
would  be  required  to  encircle  the  Z  axis  and  return  to  the 
transverse  starting  coordinates.  On  the  average  these 
extrapolations  Indicated  that  It  takes  about  100,000  first 
order  cycles  for  the  field  lines  to  encircle  the  axis. 


We  have  Investigated  the  magnetic  field  of  the 
dlsplaced-coil  configuration.  Two  characteristics,  the 
uniformity  of  field  magnitude  and  the  ratio  of  transverse 
to  axial  field  components,  have  been  analyzed  as  functions 
of  spaoe,  dlsp,  alpha,  and  location  within  the  field.  The 
pattern  of  field  lines  has  been  ascertained. 

The  qualities  of  a  magnetic  field  which  are  Important 
to  plasma  oontalnment  seem  to  be  controllable  by  adjustment 
of  the  coll  configuration.  The  field  lines  are  helices. 
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and  the  entire  field  seems  to  rotate  about  the  central 
axis . 

Note  that  the  dlsplaced-coll  configuration  does  not 
depend  on  any  particular  solenoldal  windings  as  far  as  shape 
or  size  are  concerned.  Standard  solenoids,  commercially 
available,  could  be  used  as  long  as  the  minimum  required 
field  strengths  could  be  attained.  Field  design  Is  accom¬ 
plished  solely  through  the  geometric  arrangement  of  the 
colls.  Thus  the  necessity  of  designing  and  constructing 
special  windings  Is  eliminated.  Many  different  configura¬ 
tions  could  be  experimentally  tested  cheaply  and 
conveniently . 

It  Is  mentioned  In  Chapter  I  that  stellarator-type 
revolutions  pf  the  field  lines  can  be  useful  In  cancelling 
the  drifts  of  charged  particles  In  a  toroidal  magnetic 
field.  Although  this  study  concentrates  on  an  Infinite 
straight  coll  configuration,  It  was  demonstrated  that  rota¬ 
tion  of  the  field  lines  does  take  place.  The  toroidal  coll 
configuration  from  which  the  straight  system  Is  derived 
would  also  cause  the  field  lines  to  rotate.  Hence  the 
advantages  of  the  stellarator  would  be  available  without 
the  additional  stellarator  windings. 
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III.  Plasma  In  Equilibrium 

In  this  chapter  the  steady-state  hydrodynamic  equations 
of  a  plasma  In  equilibrium  In  an  axlsymmetrlc  system  are 
derived.  The  distribution  function  for  a  system  of  like 
particles  Is  first  derived  from  Llouvllle's  equation.  The 
distribution  function  Is  applied  separately  for  Ions  and 
electrons.  The  hydrodynamic  equations  for  each  distribution 
are  found  as  moments  of  the  Boltzmann  equation.  Two  of 
Maxwell's  equations  can  be  applied  to  the  axlsymmetrlc 
system.  The  result  Is  a  system  of  eight  simultaneous,  first 
order  differential  equations  in  terms  of  eight  unknown 
variables. 

Distribution  Function 

From  Llouvllle's  theorem  (Ref  7 156)  the  most  general 
equilibrium  distribution  for  a  system  of  like  particles  Is 
a  function  only  of  the  Hamiltonian  constants  of  motion.  In 
an  axlsymmetrlc  system  two  appropriate  dynamical  constants 
are  the  Hamiltonian  and  the  canonical  angular  momentum. 

H  -  'k.mv1  Cj^<J>  (2) 

-  rmrui  -+■  c^r  A®  (3) 
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Then  a  distribution  function  which  is  consistent  with  the 
conservation  of  H  and  is 


■f1  -  £0  e  X  p  —  /3  [h 

=  -P0ex^>  —f3  [imv1  +  <^-XLmruj -_n.C|rAve] 


(4) 


where  f„  ,  and  Cl  are  arbitrary  oonstants.  In  prlnolple 

all  variables  of  the  system  of  particles  obeying  this 
distribution  can  be  found  as  moments  of  the  function.  How¬ 
ever,  it  is  convenient  to  transform  to  another  coordinate 
system  before  taking  moments. 

Let  a  coordinate  system  be  chosen  in  which  the  orthog¬ 
onal  unit  vectors  are  8^,  &2,  and  such  that 

B  =  B£.  .  (5) 

and 

§x=  ©  x  (6) 


The  new  system  is  illustrated  in  Pig.  7. 

Let  us  make  the  convention  that  the  ®  direction 
always  oolncldes  with  the  azimuthal  direction  relative  to 
the  original  cylindrical  coordinate  system.  Then  the  prod¬ 
uct  (r ck>)  in  Eq  (3)  is  uzlmuthal  velocity,  and  is  in  the 
@  direction.  The  other  two  components  of  total  particle 
velocity  can  be  called  V-^  and  Vg.  Thus 

^  v*  +•  [ruj]1  (?) 
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Fig.  7.  Cylindrical  Coordinate  System 
The  distribution  can  now  be  written 

£  =  exp  ~/3  - r_n.A©] 

■i-'lim  [v*  -4-  \£  -+-  [ru)]x]  —  rxumrujj 

If  we  add  and  subtract  the  quantity  ^mCrU)2  in  the  expo¬ 
nent,  the  terms  of  f  can  be  rearranged* 

P  =  -PQ  exp  -/#[cjJ(J>-r_n_Ae]  + 

-am[r\nja  +■  -irm[ruj  -r_n_]:LJ 


(a) 


(9) 


Certain  moments  can  now  be  conveniently  calculated. 


Particle  density  Is 


GSP/PH/69-7 


where  each  component  of  velocity  must  be  integrated  from 
minus  infinity  to  plus  infinity. 

(ID 

(Ref  8i284).  Define  the  average  azimuthal  particle  veloc¬ 
ity  as  W. 


Average  particle  kinetic  energy  is 

X  =  2  r?" jV-Pdv  =  qs  +  T  "M"  (13) 

These  last  three  equations  suggest  the  following 
conclusions i 

1.  The  density  distribution  is  not  necessarily 
constant  with  radius  in  the  azimuthal  system. 

2.  The  plasma  exhibits  rotation  as  a  solid  at  an 
angular  frequency  .CL. 

3.  Particle  energy  is  partitioned  between  the 
kinetic  energy  associated  with  random  velocities  and 
the  kinetic  energy  of  an  average  rotation.  If 

@  1/KT ,  the  distribution  behaves  like  a  Maxwellian 
distribution  with  a  drift  velocity  of  (rXL)  super¬ 
imposed. 

The  distribution  function  derived  above  describes  a 
system  of  like  particles  in  an  axial  magnetic  field.  It 
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Implicitly  accounts  for  Interactions  between  like  particles. 
The  derivation  Is  not  valid  If  different  types  of  particles 
are  considered  together  because  It  does  not  account  for 
collisions  between  unlike  particles.  However,  a  real  plasma 
composed  of  weakly  Interacting  electrons  and  ions  can  be 
described  by  applying  Eq  (9)  separately  for  each  distribu¬ 
tion  if  the  interaction  between  them  can  be  expressed.  The 
two  distributions  can  be  approximately  connected  via  an 
Inter- particle  collision  force.  Appendix  D  derives  the 
expression  for  such  a  force.  The  Interaction  force  can  be 
expressed  as  an  average  force  arising  from  collisions  be¬ 
tween  electrons  and  Ions.  The  force  Is  proportional  to  the 
product  of  the  two  local  densities  and  the  difference 
between  the  two  average  rates  of  rotation.  The  force  then 
can  be  thought  of  as  an  average  drag  between  two  distribu¬ 
tions  rotating  at  different  frequencies. 

Hydrodynamic  Equations 

This  section  outlines  the  derivation  of  the  equations 
describing  a  system  of  electrons  and  Ions  nearly  in  equi¬ 
librium.  The  relations  employed  are  the  Boltzmann 
phase-space  equation  and  Maxwell's  equations.  Mks  units  are 
used  throughout. 

Assumptions .  The  assumption  of  an  axially  directed 
magnetio  field  has  already  been  made.  Asstime  also  a  purely 
radial  electric  field.  It  was  mentioned  In  Chapter  I  that 
we -wish  to  Investigate  the  effects  of  artificially  Inserting 
Ions  and  electrons  into  the  plasma.  Assume,  therefore,  that 
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particles  may  be  inserted  anywhere  within  the  plasma  at  any 
time  rate  desired.  For  instance,  ions  and  electrons  may  be 
deposited  separately  or  together,  according  to  some  contin¬ 
uous  radial  distribution  or  perhaps  only  within  a  cylindrical 
band  of  some  /\r  in  the  system. 

Boltzmann’s  Equation  and  Moments .  The  steady-state 
Boltzmann  relation  for  a  system  of  like  particles  obeying 
the  distribution  function  of  f  is 


V  •  Vp  +  F  •  vv-P  =  b£\ 

m  a+|coll 


(14) 


where 

■f  $  (r,y) 


(15) 


F  =  F  (r,  v) 


(16) 


(Ref  I61I55).  F  is  used  here  as  the  sum  of  all  applied 
forces  which  may  be  described  as  averages.  df/9t  accounts 
for  processes  such  as  ionization  or  recombination,  charge 
exchange,  or  creation  of  particles. 

F  is  composed  of  the  following  1 

F  =  Q  [E  +  V  X  B]  +  Fs  -V  Ec ./n  (17) 


The  first  term  is  simply  the  classical  Lorentz  force  on  a 
charged  particle  of  velocity  v.  Fc  is  the  average  volu¬ 
metric  force  due  to  collisions  between  unlike  particles. 
It  is  divided  by  density,  n,  in  order  to  reduce  it  to  a 
particulate  term. 
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Fs  Is  the  average  force  per  particle  necessary  to 
accelerate  newly  created  particles  from  their  original  input 
velocities  to  the  local  average  velocity  of  like  particles 
in  the  plasma.  The  process  of  acceleration  is  assumed  to  be 
by  elastic  collisions  with  like  particles  only.  This 
assumption  is  Justified  by  the  very  small  degree  of  energy 
transfer  in  collisions  between  electrons  and  ions.  Electron- 
ion  collisions  tend  mainly  to  randomize  the  electron 
velocities!  ion  velocities  are  not  appreciably  changed  in 
such  collisions. 

Define  S  as  the  time  rate  of  creation  of  particle 
density  at  any  given  point.  Consistent  with  this  definition 

£s  =  -m  [y  -  Yo]  -S-  (18) 

where  v  is  local  velocity  and  v©  is  the  original  velocity 
of  the  new  particle.  The  minus  sign  on  the  force  conforms 
to  the  use  of  Fs  as  a  force  applied  to  the  distribution. 

Let  us  assume  henceforth  that  new  particles  are  input  at 
zero  velocity.  This  is  not  a  dangerous  assumption  in  that 
we  must  assume  that  we  may  control  the  input  velocity  any¬ 
way.  It  will  be  seen  that  v0  can  be  reinstated  at  a  finite 
value  at  a  later  point  in  the  derivation  with  little 
difficulty . 

Macroscopic  equations  of  a  plasma  are  derived  as 
moments  of  the  Boltzmann  equation.  The  first  moment  is  the 
continuity  relation. 
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V  •  [^n]  =  S  '  (19) 

The  second  moment  describes  momentum  transfer  per 
particle.  In  steady-state i 

* — » 

m  [y  •  v]  v  =  -  v-  P  +F  (20) 

n 

(Hef  l6tl6l),  where  P  Is  the  pressure  tensor. 

If  a  Maxwellian  distribution  Is  assumed  locally  for 
each  distribution,  the  pressure  tonsor  Is  a  diagonal  with 
equal  terms  (Ref  19«24). 

P  =  n  K  T  (21) 

The  most  Important  assumption  Implicit  In  the  use  of  this 
term  is  that  the  pressure  on  a  single  distribution  Is' due 
only  to  random  motion  of  particles  and  collisions  between 
like  particles. 

Maxwell’s  Sq uat 1 ons .  Under  the  assumptions  In  force 
two  of  Maxwell's  equations  supply  useful  relations. 

V  X  B  =  j  (22) 

v  -_E  =  ^4  (23) 

where  Is  charge  density  and  j  Is  current  density. 

Vector  Component  Eq uat 1 ons .  In  cylindrical  coordi¬ 
nates  the  continuity  relation  Is  simply  a  scalar  equation. 
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J_  d[rnvr]  _  ^ 
r  dr 


(24) 


where  vr  Is  the  average  radial  drift  velocity  of  particles 
out  of  the  system.  Note  that  only  partial  derivatives  with 
respeot  to  r  are  non-zero  in  an  axlsymmetrlc  system. 

The  momentum  transfer  vector  equation  supplies  two 
vector  component  equations.  The  radial  vector  component 
equation  is 


m 


[Vr  dvr  ~  Me1"]  —  ~KT  dn  ~V~  Q  ~TnMr  S 

L  dr  r  J  n  dr  r  n 


(25) 


(Ref  19i219)  where  ve  is  the  average  azimuthal  velocity  of 
particles  about  the  ax's  of  the  system.  The  azimuthal 
vector  component  equation  is 


m 


f  .  d  Vr-  V»  "|  _  P\  S 

[Vrdr  +  — J-  -^vrB-mve—  ±  — 


(26) 


Fc  appears  only  in  this  equation  because  it  is  azimuthal  in 
direction  only.  Hence  we  write  It  as  a  scalar  in  this  equa¬ 
tion.  The  sign  of  Fc  depends  on  whether  Eq  (26)  is  writter. 
for  electrons  or  ions.  Let  the  electron  distribution  be 
chosen  as  the  zero-subscripted  distribution  in  Appendix  D. 
Then  Fc  should  appear  as  a  negative  term  in  the  electron 
equation  and  as  a  positive  term  in  the  ion  equation. 

In  an  axlsymmetrlc  system  Eq  (22)  becomes 
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dB  _ 

dr  ^°Je 


The  azimuthal  current  density,  J#  ,  can  be  written  so  that 


dB 

dr 


"  [n<oVe  -  Zn,UL, 


J 


(28) 


where  e  Is  the  electronic  charge.  The  subscripts  e  and  1 
are  introduced  to  differentiate  between  the  electron  and 
Ion  distributions.  Henceforth  v  will  be  used  to  refer  to 
electron  velocities  and  u  to  refer  to  Ion  velocities. 


The  electric  field  equation  Is 

J_  d  rE  =  e_[Zn;~nt] 
r  dr  e0 


(29) 


where  we  have  rewritten  charge  density  as  the  sum  of  the 
electron  and  Ion  charge  densities. 

We  have  then  a  set  of  simultaneous  first  order  differ¬ 
ential  equations  with  r  as  the  Independent  variable.  There 
are  two  momentum  transfer  equations  and  one  continuity  equa¬ 
tion  for  each  particle  distribution,  In  addition  to  two 
Maxwell  equations  for  B  and  E,  for  a  total  of  eight  equa¬ 
tions.  The  equations  are  In  terms  of  particle  densities, 
average  radial  and  azimuthal  velocities,  radial  el«;ctrlc 
field  strength,  and  axial  magnetic  field  strength. 

These  variables  are  the  unknowns  of  interest.  We 
wish  eventually  to  study  the  radial  profiles  of  each.  This 
does  not  necessarily  mean  that  we  must  solve  each  of  the 
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eight  equations  for  one  of  these  unknowns.  Several  possible 
changes  of  variables  are  considered  In  the  next  chapter. 

Normal lzatl on.  Before  moving  on  to  the  next  chapter, 
wherein  attempts  are  made  to  solve  the  eight  equations  Just 
derived,  we  Introduce  normalization  of  the  variables. 
Normalization  eliminates  dimensional  units  and  provides  very 
convenient  scaling  o*'  our  particular  variables. 

For  instance,  let  all  velocities  be  normalized  In 
terms  of  the  mean  random  velocity  of  the  positive  Ion: 


v  = 


3KT 

mi 


(30) 


where  1^  Is  temperature  in  degrees  Kelvin.  The  object  Is  to 
replace  all  velocities  with  corresponding  primed,  dimension¬ 
less  numbers’.  For  convenience  define 


(31) 


Let  time  be  normalized  In  terms  of  the  cyclotron  period  of 
the  Ion  near  the  axis. 


-t  -t'T 


(32) 


T  = 


m, 

ZeBo 


(33) 


where  B0  Is  the  axial  magnetic  field  strength.  It  was 
assumed  earlier  that  B0  could  be  specified. 

Earlier  it  was  mentioned  that  one  of  the  goals  of  this 
study  is  to  determine  theoretically  the  feasibility  of 
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attaining  particle  number  densities  on  the  order  of 
1020  m“^.  Let  densities  and  the  density  source  terms  be 
normalized  to 

20  -3 

Tl  =  10  m  04) 

Consistent  with  the  above  definitions  are  the  follow- 
1 ng  normalizations! 


Radius i 

r  =  r'VT 

(35) 

Magnetic  field i 

B-  B'B0 

(36) 

Rate  of  particle  ere at 1 on « 

c  S’n 

(37) 

Angular  frequency i 

% 

it 

3 

(38) 

Electric  field  i 

(39) 

Collision  forcei  C 

■c 

r-<  "V 

:  “  *c  'Y  mi 

(40) 

When  these  definitions  are  substituted  into  the  eight 
equations  all  dimensions  cancel  out.  Let  us  then  drop  all 
primes,  remembering  from  the  above  definitions  how  to 
recover  the  real  variables.  Table  I  is  a  compilation  of  all 
eight  normalized  equations,  plus  Fc ,  the  normalized 
collision  force  from  Appendix  D,  after  the  primes  are 
dropped. 


# 
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Table  I 


1  1 C fin*oz)  -  Sc 
*  J~7L 

i  cHoJli  =  s'/ 
a  cfiL 


(41) 

(42) 


me(iTrJirr  -u^J  _  -  rru  Te  _£  fir*  g)  5e 

\  *  J  3  Ts  ncJZ  z  7?7 

»7?,7#r  JUr-Ka)  -  -J7U-  dri:  +  mi(E^oS)  -  S; 

\  J X.  n)  3  n;  d  n.  W7 


(43) 

(44) 


iTr^/^  +  Jg]  -  \Tr  B  -  nit.  Vi  5J:  - 

\Jn  nj  Z  /?e 

2 CrfJ  Mt,  ±  Ub\  —  —  Mi  2Cr  B  -  IT)  i  lA-oSi  f. 

'  U  ■*■  * )  n7  n7 


d£  x  3JLTi  (nK\ra  -  2 /?/ zce) 
ZeB0  K  ' 


±  d(n.E\  =  mi  n 
cl  !L  F^BT 


(45) 

(46) 


(4?) 

(48) 


(49) 


(50) 
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Summary 

This  chapter  is  a  derivation  of  a  set  of  time- 
independent  differential  equations  which  describe  the 
macroscopic  variables  of  a  system  of  ions  and  electrons 
near  equilibrium.  It  is  useful  to  summarize  the  assumptions 
upon  which  the  derivation  is  based i 

1.  The  system  is  axisymmetrlc. 

2.  Both  ions  and  electrons  obey  Boltzmann-like 
distributions  and  a  real  plasma  can  be  accurately 
described  by  applying  their  respective  distribution 
functions  separately. 

3.  The  magnetic  field  is  axially  directed  and 
its  strength  on  axis  can  be  specified. 

4.  The  electric  field  is  radial  in  direction. 

5.  Ions  and  electrons  can  be  Inserted  into  the 
plasma  arbitrarily  at  any  time  rate  desired.  Decel¬ 
erating  drag  forces  on  the  two  separate  distributions 
arise  from  collisions  between  like  particles  only,  and 
may  be  treated  as  average  forces. 

6.  The  drag  force  between  electron  and  ion 
distributions  may  be  treated  as  an  average  volumetric 
force . 

7.  Pressure  on  a  single  distribution  is  due  only 
to  random  collisions  between  like  particles  in  the 
distribution. 
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8,  The  electron  and  Ion  temperatures  are  known 
and  are  uniform  across  the  plasma. 
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IV.  Solving  the  Differential  Equations 

Introduction 

This  chapter  describes  In  detail  attempts  to  develop 
procedures  for  the  numerical  Integration  of  the  system  of 
differential  equations  derived  in  Chapter  III.  The  objec¬ 
tive  of  these  attempts  was  to  select  In  terms  of  accuracy 
and  efficiency  the  best  methods  for  numerical  solution  of 
this  particular  system  of  equations.  Involved  are  the 
choice  of  variables  sought  as  solutions,  the  adaptation  of 
standard  numerical  Integration  techniques  to  the  specific 
problem,  and  the  writing  of  computer  programs  to  perform 
calculations  most  accurately  and  efficiently. 

First  a'  set  of  variables  Is  selected  and  the  equations 
In  Table  I  are  solved  for  the  first  derivative  of  one  of 
the  variables.  Equations  for  the  Initial  conditions  on 
the  variables  are  derived  from  the  equations  written  on 
axis.  Sample  calculations  are  made  of  non-zero  initial 
conditions.  These  representative  values  for  the  Initial 
conditions  serve  to  point  out  potential  sources  of  error 
In  the  computation  of  the  derivatives.  Various  changes  of 
variables  are  studied  In  attempting  to  eliminate  the 
sources  of  error.  The  chapter  Is  concluded  with  an  anal¬ 
ysis  of  the  application  of  numerical  integration  techniques 
to  the  solution  of  the  differential  equations. 
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Choloe  of  Variables 

In  Chapter  III  mention  Is  made  of  which  of  the  vari¬ 
ables  Involved  in  the  system  of  equations  are  of  Interest, 
but  equations  are  not  selected  to  be  solved  for  specific 
variables.  In  this  section  alternative  methods  of  solving 
the  equations  are  discussed. 

The  following  macroscopic  quantities,  as  functions  of 
radius,  are  the  objects  of  this  study  ■  ne ,  nt ,  vr,  ve  ,  u0, 
B  and  E.  (The  equations  In  Table  I  are  In  terms  of  the 
corresponding  normalized  variables,  but  the  conversions 
back  to  the  real  quantities  are  simple.)  A  straightforward 
approach  Is  to  solve  each  equation  for  the  first  derivative 
with  radius  of  each  of  these  unknowns,  but  other  arrange¬ 
ments  of  the  equations  are  possible  and  may  be  more  useful. 
This  chapter  is  primarily  concerned  with  evolving  a  set  of 
variables  and  the  solutions  for  their  first  derivatives 
which  are  most  amenable  to  numerical  Integration,  keeping 
in  mind  the  limitations  Imposed  by  the  available  numerical 
techniques  and  computer  technology. 

Eq  (47)  for  magnetic  field  In  Table  I  is  apparently 
already  In  Its  most  useful  form.  There  are  no  obvious 
rearrangements  which  might  improve  it. 

Consider  Eq  (45).  Move  the  second  term  on  the  left 
hand  side  to  the  right  hand  side  and  divide  through  by 
mevr.  The  result  is  a  direct  expression  for  dvQ/dr. 
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d  sfe  =  _  _v»  _ 

rSe  ■+■  m; 

-  -M 

dr  r 

However,  note  that 

vr  me 

L2 

neVrJ 

d  Ve  Va  =  _l_  d[rVal 

dr  r  r  d  r 


Thus,  rearranging  Eq  (45) 


r  S-  +  m;  B  _  FL 

Vr  ne  me.  |_  2  Mr 


(5D 


(52) 


(53) 


The  Independent  variable,  r.  Is  always  known  exactly  In  the 
step-wise  numerical  Integration  techniques  usually  employed 
with  first  order  differential  equations.  Therefore,  there 
should  be  no  trouble  retrieving  ve  from  the  solution. 
Similarly  Eq  (46)  becomes 


d  ruu  =  -  ruL»  Si  -  B  -  Fc. 

dr  u.r  n;  [_  rv.  u.r. 


(54) 


Divide  Eq  (41)  by  n@vr  and  expand  the  derivative. 


I  dne 

He  dr 


_J _  d  Vr-  = 

Vr  dr 


He  Vr 


Note  that  this  equation  could  lie  solved  for  the  first 
derivative  of  either  nQ  or  vr<  For  the  moment  use  Eq  (55) 
to  eliminate  dvp/dr  from  Eq  (43).  Then  a  solution  for 
dne/dr  Is 
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Vo  4-  Vr 

1*  Vr  —  Se. 

-  m. 

[E-t  Bf] 

=  ne 

.  r 

L  r  Oa. 

2.  me 

L  el 

f  )  Te.  m  i  -  Vrx 

[3  Ti  me 


Similarly  for  the  ion  equation 

r. .  X  . .  r. . 


1 

=.  nt  L 


U-r  ttr  —  3-Si  1  -v-  E  -v-  iLefi 
L  r  n.  J  L _ 

",  i  -  Ul^"’ 

.  3 


Assuming  that  ng  is  a  wise  choice  of  variable,  Eq  (55) 

could  be  solved  for  dv  /dr  and  the  result  integrated  for  v 

v  r 

directly.  However,  if  rnevr  is  chosen  as  a  variable, 
d(rnevr)/dr  is  rSe ,  a  quantity  which  has  been  assumed  to  be 
exactly  specified.  Therefore  the  differential  of  rnevr  is 
known  exactly,  and  for  this  reason  it  seems  a  wise  choice 
of  variable.  Hence, 

rna  Vr 

vr  -  rne  (58) 


rn.  u.r 
°-r  =  rn. 

Eq  (48)  can  be  left  almost  as  is, 

d[rE]  =.  mi  r  [iErv,  —  r\eJ 

j -7  r~>2_  L  J 


dr  ZeA  L 

The  equations  in  Table  I  of  Chapter  III  have  been 
rearranged  as  solutions  for  the  first  derivatives  of  the 
following  variables!  rnevr»  ra^,  rv4  ,  ru,,  B,  and  rE. 
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The  next  logical  step  is  to  find,  the  Initial  conditions  on 
all  variables  .so  that  the  differential  equations  can  be 
solved. 

Initial  Conditions 

If  the  numerical  Integration  of  the  system  of  differ¬ 
ential  equations  Is  to  be  started  from  zero  radius.  It  is 
necessary  to  have  a  means  of  computing  the  numeric  values 
of  all  variables  on  axis.  Some  of  the  variables  are  zero 
by  symmetry i  those  others  which  are  not  arbitrarily  con¬ 
trolled  must  be  analytically  derived.  Then  standard 
numerical  methods  of  solving  first  order  simultaneous 
differential  equations  can  be  used  to  carry  on  the  integra¬ 
tion  away  from  the  axis. 

Some  of  the  initial  conditions  are  immediately  evident 
from  symmetry.  Recall  that  the  radial  and  azimuthal  veloc¬ 
ity  terms  are  averages,  or  drift  velocities,  referenced  to 
the  axis.  Then  on  axis  all  the  velocities  are  zero.  The 
electric  field  must  be  zero  on  axis  by  symmetry  also. 

Although  the  azimuthal  velocities  at  r=0  must  be  zero, 
the  angular  frequencies  of  rotation,  COe  and  CU;  ,  can  be 
finite.  Indeed  In  Chapter  III  it  was  suggested  that  such 
frequencies  are  constant  across  the  plasma. 

Since  in  Chapter  III  the  magnetic  field  Is  normalized 
to  its  axial  magnitude,  B  is  unity  on  axis.  Recall  the 
assumption  that  the  real  value  of  B  on  axis.  Bo  ,  can  be 
arbitrarily  Imposed.  It  was  also  assumed  that  the  two 
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distribution  temperatures,  Tg  end  ,  are  known  and  are 
uniform  across  the  plasma. 

The  electron  and  ion  number  densities,  n  and  n  ,  are 

6  1 

the  remaining  initial  conditions.  Let  us  seek  algebraic 
solutions  for  and  n.  by  solving  simultaneously  the 
differential  equations  as  they  can  be  written  or.  axis. 
First  multiply  Eq  (55)  by  r. 

r  -4-  r  dne.  -+•  jr_  dvr  -  _r_  ik 
r  ne  dr  vr  dr  vr  rie 


(6i) 


Take  the  limit  of  vr/r  as  r  goes  to  zero 


JLurrv  Vf  =  dvrxdr  =  dvr 

r— >  o  r  dr/dr  dr 


(62) 


(Ref  14 « 262 )_ .  Since  dne/nedr  must  be  finite  on  axis,  the 
second  term  in  Eq  (61)  is  zero  on  axis.  Then 

0_  JVr_  =  -Sei 

r  ne 


(63) 


Reduce  Eq  (56)  to  an  algebraic  equation  with  no  vr  or  vc 
terms  by  use  of  Eqs  (6l),  (62),  (63),  and 

Ve 


UUe  = 


(64) 


The  result 

_a. 

is 

LSe' 

1 

1 

=  UJe.  -  ml 

Ul)»  -4-  C 

M- 

.  He  _ 

£me 

r  _ 

(65) 


where  B  has  been  set  equal  to  one.  Similarly  for  the  ion 


distribution 
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3 

4- 


1. 

U>i 


u);  -t- 


r 


(66) 


At  this  point  It  Is  convenient  to  make  a  change  of 
variable  for  the  purposes  of  the  solution  for  initial  condi¬ 
tions.  Define  as  the  difference  between  densities 

f0  =  2n;-ne  (67) 


Also  define  an  average  density 

n  =  [Z  n,  4-  nej/a 

As  In  Eq  (62)  when  r  goes  to  zero, 

JL™  _E  - 

o  r  dr 

Therefore  Eq  (60)  becomes  on  axis 

Oi  E  =  rrv,  /° 
“r  ZeB0x/ 

Define  for  subsequent  use 


(68) 


(69) 


(70) 


r  _  nn, 

*  ZlB?  <: 

Turn  to  Eqs  (53)  and  (54)  to  see  If  expressions  for 
{U,.  and  LO\  can  be  derived.  Using  the  definition  of  UJe  , 
Eq  (64),  and  the  Initial  value  of  B,  Eq  (53)  becomes 
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HuJe  =  -2.UJe  +  _DQl 

me 

Assume  In  the  normalized  expression  for  Fc  In  Table  I  that 
the  first  power  of  the  difference  (Cu>e  -  UJj  )  Is  much 
greater  than  any  of  Its  higher  powers.  This  assumption 

will  be  verified  when  GJe  and  UJ',  are  solved  for.  Then 

/ 

approximately 

Fc  =  Cf  zEn.IHe  [uie  -U)»]  (73) 


J_  -  Fc 
Z  n<s.vr 


(72) 


includes  all  the  constants  in  the  first  term  of  the 
series  of  the  normalized  expression  for  Fc.  A  glance  at 
Appendix  D  shows  that  in  a  strict  sense  an  unknown,  ln(ng), 
is  thus  Included  in  Cf.  However,  it  happens  that  the 
normalization  of  densities  causes  ln(ng)  on  axis  to  pe  very 
near  zero.  When  the  actual  method  of  computing  the  axial 
ne  and  n^  is  described,  it  will  be  shown  how  the  value  of 
lnfrig)  Is  corrected  in  the  coefficient  Cf. 

Using  the  change  of  variable  previously  defined. 


Zn,r\£  - 


(74) 


Assume  for  the  moment  that  the  relative  difference  between 
Zn^  and  rig  is  so  small  that 

"  C7; 

This  approximation  will  be  validated.  Therefore 
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(76) 


(77) 


Similarly  from  Eq  (5*0  on  axis  is 


~J-  _u  ^ 

rrVi  C-P 

1  1 
.ZS-,  's: 

na 

8  JDs.  4-  4-Cf 
m, 

i  4-  1 

.  Se.  rn  i  <8 1  . 

n2- 

The  assumption  made  in  Eq  (75)  also  allows  the 
following i 

J _  ^  I 

He1  [na  -/^n]  (79) 


ZniHe  — 


Then  using 

Eq 

(63) 

s  olve  Eq 

(53) 

for  dJ, 

a. 

4- 

r  1 

I  1 

uJ/o  =  Z. 

C-f 

S=.J 

8  DDfi.  4-  4- C-f  r_L_  4-  me-  _L_1  nx 

rr\t  ISe.  m,  5.  - 


Expand  this  result  in  a  binomial  series 


/£  3  /■==>  a 

I  =14-  n  4-4  n  4  •  •  • 


°x  [ 

n3" 

i 

(80) 

[ 

n  J 

Invoking 

Eq  (75)  again. 

1  ~  i 

r\tx  n2- 

\  4- 

n 

(81) 

Similarly  for 
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Degree  of  Term 


T Sable  II 

Coefficient 
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I 


(82) 


2  2 

These  approximations  can  "be  used  to  replace  n_  and  n.  in 

e  1 

Eqs  (65)  and  (66).  Eq  (70)  gives  E/r  In  terms  of  and  a 

constant.  Thus  Eqs  (65)  and  (66)  can  both  be  solved  for 
In  terms  of  n,  COe  ,  CO;  ,  Se ,  and  S1#  Eq  (65)  Is  solved 
for  ^  below  as  an  example . 


- f-Sq-.  1  4-  UJeX 

^  LnJ _ 

_3_  ■Se.2,  4-  m  .  £ 

4-  n3  Z  me. 


m.  uJe 
,Z  me. 


(83) 


If  the  solutions  for  are  equated,  a  ninth  degree  poly¬ 
nomial  In  n  can  be  written.  Table  II  Is  a  list  of  the 
coefficients  of  each  power  of  n. 

There  are  methods  of  solving  for  the  roots  of  poly¬ 
nomials  of  any  degree  (Ref  IO1I6).  However,  it  Is  not 
necessary  to  solve  for  all  nine  roots  of  the  polynomial. 

We  are  interested  only  in  those  positive  roots  which 
correspond  to  an  unnormalized  density  on  the  order  of 
102Om-3.  It  can  be  shown  that  these  roots  can  be  found 
directly. 

Approximate  Solution.  Instead  of  trying  to  solve  the 
entire  ninth  order  polynomial  for  all  nine  roots  of  n  let 
us  seek  an  approximate  solution  for  the  root,  or  roots,  of 
interest.  Note  that  the  coefficients  in  Table  II  are  in 
terms  of  Cf ,  £  ,  Se ,  and  Si.  Cf  and  <5  are  the  normalized 
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coefficients  of  the  equations  for  Fc  and  d(rE)/dr, 
respectively.  Se  and  Si  are  the  normalized  source  terms  on 
axis.  It  is  very  Interesting  to  invent  a  numerical  example 
which  demonstrates  how  the  solution  for  n  depends  on  these 
s  ource  te  rms . 

Suppose  that  the  unnormalized  source  terms  are  on  the 
order  of  101^m“-^/sec.  That  is,  if  we  begin  with  no  density 
at  all  and  assume  no  diffusion  away  from  the  region  of  the 
axis,  in  ten  seconds  the  number  densities  of  both  electrons 
and  ions  would  be  about  102®m“3.  Suppose  also  that  Se  and 
Si  differ  by  one  percent  in  magnitude .  At  this  point  we 
will  not  specify  which  is  larger.  Let  B0  equal  one  kilo- 
gauss.  If  now  the  source  terms  are  normalised,  and  put 
into  the  expression  for  the  coefficients  in  Table  II,  the 
non-zero  coefficients  assume  the  following  approximate 
value  s i 


Degree  of  Coefficient 
9 
7 
6 
5 
4 
3 
2 
0 


Value  of  Coefficient 
-1012 
1012 
10”  9 
106 
10”3 
10-13 

lO”1^  . 

10-37 
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Note  that  the  ninth  degree  coefficient  is  the  only 
negative  coefficient.  An  examination  of  Table  II  shows 
that  this  coefficient  must  always  be  negative.  The  other 
coefficients  are  all  positive.  Since  only  positive  roots 
of  the  polynomial  interest  us,  the  magnitude  of  the  ninth 
degree  term  of  the  polynomial  must  always  equal  the  magni¬ 
tude  of  the  sum  of  the  other  nine  terms  if  the  polynomial 
is  to  sum  to  zero.  For  the  set  of  values  Just  computed  it 
is  evident  that  the  seventh  and  ninth  degree  terms  are  by 
far  the  largest  terms  in  the  polynomial.  To  illustrate, 
suppose  we  approximate  the  polynomial  as  follows  i 

-10V+  io'V-o  <84) 

The  solution  is  obviously  n  equal  to  about  ons.  Hence  for 
this  root  the  other  terms  of  the  polynomial  are  indeed 
small  by  comparison  and  can  be  neglected.  This  root  trans¬ 
lates  to  an  unnormalized  average  density  of  about  102®m"^, 
which  is  about  the  size  of  number  densities  of  interest. 
Therefore  the  source  terms  chosen  for  this  example  at  least 
fall  in  the  neighborhood  of  Interesting  oases. 

For  specified  conditions  very  different  from  those  of 
the  above  example,  the  approximations  made  should  be 
checked  for  validity.  Nevertheless,  if  the  approximations 
do  hold  nearly  as  well,  an  approximate  expression  for  n2  is 
then  available  from  the  ratio  of  the  seventh  to  the  ninth 
degree  coefficients  in  Table  II.  It  is  interesting  to 
examine  this  ratio,  using  only  the  first  term  of  the 
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seventh  degree  coefficient  since  it  is  much  larger  than  the 
other  term. 


8  n  + 

nrv.  me 

L 

Zme.  m> 

C*  Se 

r 1  m 

a 

LZS;  SeJ 

(85) 


Neglect  terms  of  order  mg/m^  compared  to  one. 


n 


o. 


8Se. 

zcp  r^Se-zs-,  ]x 
zs,  J 


(86) 


Note  that  n  is  inversely  proportional  to  the  relative 
difference  between  Se  and  ZS^  and  proportional  to  the 
square  root  of  the  source  terms. 

It  Is  to  be  stressed  that  the  above  approximate 
expression  for  n~  may  not  be  valid  for  all  Imposed  initial 
conditions.  It  will  be  used  only  as  a  means  to  estimate 
the  root  of  interest.  The  complete  polynomial,  without 
approximations,  should  be  computed  to  insure  that  it  is 
near  zero.  If  It  is,  the  solution  for  n  can  be  refined  by 
Iterative  correction  schemes.  The  Newt on-Raphs bn  method  of 
finding  a  root  of  a  polynomial  is  quite  useful  and  quick 
when  a  good  estimate  of  the  root  is  available  (Ref  10:630). 
It  is  well  suited  to  iterative  computer  techniques. 

The  Cf  used  in  Eq  (86)  to  compute  the  first  trial 
value  for  n  was  computed  with.ln(n0)  equal  to  zero.  Then 
in  every  iteration  of  the  Newt  on-Raphs  on  method  ln(n)  should 
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tie  used  where  ln(ne)  appears  in  Cf.  Ln(ne)  is  not  one  of 
the  larger  terms  in  ln(-A-).  (See  Appendix  D.  )  Then  the 
very  small  difference  between  lnfrig)  and  ln(n)  cannot  have 
a  great  effect  on  the  accurate  calculation  of  Cf.  Hence 
rapid  convergence  of  the  iterative  solution  for  the  root  of 
interest  is  realized. 

Origins  of  Error 

It  is  prudent  at  this  point  to  look  at  the  physics  of 
the  situation  and  search  for  any  arithmetic  operations  in 
the  eight  equations  which  might  introduce  error.  Since  the 
intent  is  to  use  numerical  integration  techniques  to  solve 
the  equations  in  Table  I,  it  is  imperative  that  the  differ¬ 
entials  chosen  in  this  chapter  be  computed  as  accurately  as 
possible.  Therefore  it  is  interesting  to  study  the  offeots 
of  the  error  associated  with  each  term  of  the  differentials 
and  to  investigate  how  such  errors  might  be  induced. 

In  the  electric  field  equation  the  derivative  depends 
on  the  difference  between  electron  and  ion  charge  densities  i 
both  are  integrated  variables  and  thus  subject  to  error. 

Any  computer  language,  has  only  a  finite  number  of  signifi¬ 
cant  digits  it  can  assign  to  any  single  variable.  m 
Fortran  IV,  for  example,  the  maximum  number  of  digits  avail¬ 
able  to  any  real  variable  is  sixteen.  (Ref  9x6)  If  Zn^ 
and  are  nearly  equal  on  the  average,  their  difference 
may  not  appear  except  in  the  last  few  significant  figures 
of  either.  The  accuracy  of  rE  can  thus  be  severely 
limited.  It  is  expected,  in  fact,  that  the  net  charge 


GS P/PH/6 9- 7 


densities  in  a  plasma  in  equilibrium  is  nearly  zero.  The 
only  way  of  determining  if  it  is  feasible  to  calculate  rE 
by  Eq  (60)  is  to  make  some  representative  calculations  of 
the  difference  between  Zn^  and  ne . 

Recall  the  assumption  that  the  axial  magnetic  field  is 
to  be  imposed.  The  electric  field  can  arise  only  from  the 
relative  motion  between  the  two  oppositely  charged  distri¬ 
butions.  The  average  azimuthal  rotation  of  a  charged 
distribution  results  in  a  radial  force  on  the  distribution 
in  the  presence  of  an  axial  magnetic  field.  Because  the 
electrons  and  ions  are  oppositely  charged,  the  radial 
forces  thus  generated  tend  to  drive  the  two  distributions 
in  opposite  radial  directions.  A  radial  electric  field 
appears  to  prohibit  this  separation.  Thus  a  cancellation 
of  force  terms  is  Involved  in  the  terms  (E  +  vffB)  and 
(E  +  UpB).  The  potential  for  a  prohibiting  loss  of  accu¬ 
racy  in  either  or  both  of  these  Lorentz  force  terms  exists 
if  E  is  nearly  equal  in  magnitude  to  either  v$  B  or  ueB. 

Three  terms  in  the  eight  differential  equations  have 
been  identified  as  potential  sources  of  error.  It  is 
expected  that  the  calculation  of  these  terms  might  require 
computing  a  relatively  small  difference  between  two  nearly 
equal  numbers.  The  loss  of  significance  amounts  to  a  loss 
of  accuracy  in  these  differences  in  computer  calculations. 

The  severity  of  accuracy  loss  in  the  sensitive  terms 
cannot  really  be  Judged  until  accurate  numeric  values  for 
the  variables  Involved  in  each  term  are  available .  Nor  can 
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the  effects  of  the  error  in  these  terms  on  the  accuracy  of 
the  differentials  be  determined  until  the  relative  magni¬ 
tudes  of  all  the  terms  in  each  equation  can  be  estimated. 

Numerical  Example .  To  ascertain  whether  the  above 
problems  can  arise  it  is  again  useful  to  go  through  a 
numerical  example.  Suppose  that  a  cylindrical  plasma 
sustains  a  voltage  as  high  as  10,000  volts  across  a  radius 
on  the  order  of  10  cm.  Suppose  also  that  the  electric 
field  strength  is  approximately  linear  with  radius.  Then 
at  5  cm  from  the  axis  the  electric  field  strength  is  about 
1000  volts/cm.  This  is  quite  a  high  electric  field  for  an 
ionized  gas  to  sustain. 

From  Eq  (29)  the  difference  between  the  unnormalized 
electron  and  ion  densities  corresponding  to  this  electric 
field  is  about  2.2  x  lO^m”^.  If  we  asBuir  that  the  mean 
particle  density  is  as  high  as  1020m”3,  the  relative  dif¬ 
ference  between  Zn^  and  rig  is  about  10“-5.  For  smaller 
electric  fields  this  difference  is  even  smaller. 

Note  that  by  the  equations  chosen  as  solutions  for 
the  derivatives  of  rig  and  n^  the  densities  are  solved  for 
separately.  It  then  becomes  necessary  to  compute  their 
difference  in  the  electric  field  equation.  So  in  computer 
calculations,  where  the  densities  can  only  be  known  to  a 
finite  number  of  significant  places,  their  difference  is 
known  accurately  to  fewer  significant  places  than  either 
density  beoause  the  first  few  significant  places  are  can¬ 
celled  out.  In  the  example  above  five  significant  figures 
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are  lost  every  time  the  difference  Is  calculated.  For 
smaller  electric  fields  or  higher  average  densities  the 
relative  differences  between  electron  and  Ion  densities  Is 
even  smaller. 

It  Is  possible  from  the  Initial  condition  equations 
derived  earlier  to  find  an  analytic  expression  for  the 
term  (E  +  VgB)  on  axis.  Immediately,  since  B  on  axis  Is 
one,  the  term  Is  simply  (E  +  r  u>&  ) .  Let  us  find  an  expres¬ 
sion  for  this  term  on  axis  and  find  Its  size  relative  to 
v0 .  The  conclusions  can  be  qualitatively  extended  for  the 
rest  of  the  plasma  and  to  the  equivalent  term  In  the  lor. 
equation. 

By  Eqs  (70)  and  (71),  In  terms  of  normalized  variables, 

El  '  V/?  (8?) 


Then  the  term  of  interest  relative  to  v^  Is 

£,  4-  ruje  =.  f z  me  ruie-  -  ruJc  +  ru)e.1  Zr  uue. 
ruje  L  no,  V 


-  Z  me  uJe- 

m'i 


Recall  that  UJt  has  been  normalized  to  the  Ion  cyclotron 
frequency,  and  that  It  is  used  here  as  an  average  frequency 
of  rotation  for  the  entire  electron  distribution.  Hence 
the  above  Is  less  than  one.  Thus,  since  me  Is  much 

less  than  m^,  E  and  r  CUg  are  nearly  equal  in  magnitude,  and 
there  Is  a  cancellation  of  the  first  significant  figures  In 
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taking  their  difference.  Similar  results  hold  for  the 
term  (E  +  u^B). 

The  significance  of  the  above  examples  and  arguments 
is  that  there  can  easily  be  a  severe  loss  of  accuracy  in 
computing  the  first  derivatives  of  ne  and  n^ .  Recall  that 
significant  figures  were  lost  in  the  electric  field  equa¬ 
tion.  Now  we  find  that  additional  significant  figures  oan 
be  lost  in  the  Lorentz  force  terms.  Thus  ne  and  cannot 
be  accurately  computed  to  even  the  limited  number  of  places 
available  in  the  computer. 

It  seems  then  that  some  means  should  be  sought  to 
solve  for  the  difference  between  Zn^  and  rig  directly  rather 
than  subtracting  them  every  time  the  derivative  of  rE  is 
computed.  That  is,  f?  should  be  chosen  as  one  of  the  vari¬ 
ables  of  the  system  such  that  it  can  be  sr  ed  for  dltectly. 
The  next  section  discusses  such  a  change  of  variable. 

Change  of  Variable 

It  is  suggested  in  the  last  section  that  Instead  of 
integrating  for  Uq  and  n^  it  would  be  wise  to  solve  direct¬ 
ly  for  the  difference  between  them.  It  was  shown  that  a 
severe  loss  of  accuracy  could  ensue  from  subtracting  the 
two  densities.  To  solve  for  the  first  derivative  of 

d  P  «  2  drv.  _  d  r\«- 

.  dr  dr  dr  (®9) 

By  definition  of  n  we  could  choose  it  as  the  other  density 
variable . 
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dn  —  J_  Z.  dn; 
dr  H  dr 

Thus,  In  principle,  we  need  only  add  or  subtract  Eqs  (56) 
and  (57)  to  maice  the  change  of  variables.  However,  eote 
that  there  Is  no  cancellation  of  any  equal  terms  accom¬ 
plished  by  the  subtraction  of  dne/dr  from  Zdnj/dr.  Thus 
nothing  has  been  done  to  eliminate  the  large  common  factors 
shared  by  Zn.^  and  ne . 

It  has  been  shown  then  why  an  obvious  change  of  vari¬ 
ables  does  rot  remedy  the  basic  problem  with  this  set  of 
equations.  The  problem  Is  the  loss  of  accuracy  Incurred  by 
subtracting  with  the  computer  nearly  equal  terms.  The  next 
section  explains  how  this  type  of  error  makes  Impossible 
the  numerical  Integration  of  differential  equations. 

Numerical  Integration 

Many  approximate  methods  of  numerically  integrating 
the  general  first  order  differential  equation  have  been  In¬ 
vented.  The  oldest  and  simplest  are  Euler’s  method  and  Its 
modification  (Ref  15»310).  They  both  essentially  make 
linear  extrapolations  from  one  known  value  of  a  function  to 
a  new  value  by  multiplying  the  derivative  of  the  function 
by  an  Increment  of  the  Independent  variable.  More  elegant 
methods  such  as  Approximating  Polynomials  (Ref  15i358),  and 
Milne  (Ref  15«353)  employ  the  extrapolation  method  more 
cleverly,  but  the  use  of  them  can  also  be  limited  by  this 
approximation. 
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The  Milne  method  Is  an  example  of  the  general 
predictor-corrector  scheme.  The  two  baslo  equations  It  uses 
are 

Predictor,  y^,  “  >^.3  +  ^  ^Xn-a"  Yn-I  *  (91) 


Corrector,  X\+l  Xn-<  Vn-I  (92) 

where  yn+1  Is  the  value  of  the  dependent  variable  being 
sought.  The  subscripts  denote  the  number  of  steps  by  which 
the  Independent  variable  has  been  lncrenented.  The  primed 
quantities  are  the  first  derivatives  of  the  dependent  vari¬ 
able.  The  predictor  equation  provides  an  estimate  of  yn+1 
at  the  nert  Incremented  value  of  the  Independent  variable 
based  on  the'  behavior  of  the  function  at  previous  Incre¬ 
ments.  The  corrector  equation  calculates  yn+1  over  and 
over  until  successive  answers  agree  to  within  some  accepted 
error. 

Note  that  In  the  predictor  equation  the  values  of  y  at 
four  previous  Increments  are  required.  Some  starting  solu¬ 
tion  Is  usually  required  to  calculate  the  first  three  values 
of  every  dependent  variable,  and,  of  course,  the  Initial 
conditions  for  all  variables  must  be  known.  Runge-Kutta  Is 
a  very  popular  method  of  obtaining  the  three  starting 
values  from  a  knowledge  of  the  Initial  conditions;  Runge- 
Kutta  Is  essentially  derived  from  a  truncation  of  the 
Newton  series  for  forward  ?  ^erpolatlon  (Ref  15*56),  and 
therefore  Imposes  the  approximation  of  linearity. 
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The  Runge-Kutta  method  was  applied  to  the  eight, 
differential  equations  for  several  sets  of  imposed  condi¬ 
tions  to  obtain  the  first  three  solutions  for  all  the 
variables  selected  earlier,  and  the  Kilne  method  was  applied 
at  the  fourth  step  as  explained  above.  Initial  conditions 
were  computed  as  described  earlier  in  the  chapter.  As  the 
corrector  equation  was  applied  repeatedly  at  the  fourth 
step,  the  values  of  the  variables  did  not  converge,  but 
rather  increased  or  decreased  beyond  reasonable  limits. 
Within  a  few  iterations  the  two  densities  decreased  by 
several  orders  of  magnitude,  and  E  and  the  four  velocities 
increased  in  a  like  manner. 

Apparently  the  sensitive  differences  in  the  density 
equations,  and  the  electric  field  equation  were  computed 
erroneously  due  to  small  errors  in  the  terms  Involved  in 
the  differences.  The  derivatives  of  these  terms  were  then 
inaccurately  calculated  since  the  differences  occurred  in 
the  last  significant  figures  of  the  variables.  The 
repeated  application  of  the  corrector  equation  served  to 
amplify  the  errors  rather  than  improve  the  predicted  values 
of  the  variables. 

Series  Expansions .  Suppose  the  exact  Initial  values 
of  all  non-zero  variables  were  available.  It  might  then  be 
possible  to  expand  all  variables  about  the  axis  in  alge¬ 
braic  series  in  terms  of  r.  If  the  variables  could  be 
expanded  to  enough  terms  so  that  very  aoourate  values  of 
each  a  short  distance  off  axis  were  available,  the  three 
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starting  solutions  could  be  had  by  computing  each  series  at 
three  different  radii. 

The  Taylor  Series  method  (Ref  15*328)  is  often  used  to 
begin  from  initial  conditions  the  numerical  solution  of  any 
order  differential  equation.  It  is  based  on  the  expansion 
of  the  dependent  variables  into  infinite  series  in  powers 
of  the  Independent  variable. 

y;  00  ■  y,  (Xo)  +  y,'(x0)[x- xj+  . .  <93) 

where  y^(xQ)  is  the  initial  condition  on  y^  at  xQ.  The 
primary  limitation  on  this  method  is  usually  the  excessive 
labor  involved  in  solving  algebraically  for  the  higher  order 
derivations.  However,  if  (x-xQ)  is  kept  small,  a  variable 
might  be  represented  accurately  by  only  two  or  three  terms 
of  a  Taylor  series. 

The  usual  procedure  of  computing  the  coefficients  of 
tne  (x-x0)  terms  is  to  substitute  the  truncated  series  for 
each  variable  into  the  set  of  differential  equations  and 
solve  simultaneously  for  the  coefficients.  Each  differ¬ 
ential  equation  reduces  to  an  algebraic  equation  which  can 
be  subdivided  into  equations  of  the  coefficients  of  like 
powers  (x-xQ).  There  are  always  as  many  coefficient  equa¬ 
tions  as  coefficients  if  there  is  a  differential  equation 
for  each  variable.  In  principle  the  set  of  simultaneous 
coefficient  equations  oan  always  be  solved,  often 
numerically,  for  the  numerical  values  of  the  coefficients. 
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Assume  then  that  It  is  possible  to  expand  each  of  the 
variables,  using  n  and  ^  as  density  variables,  into  series 
whloh  represent  the  variables  a  short  distance  from  the 
axis  to  greater  accuracy  than  the  digital  computer  is 
capable  of  carrying.  Each  variable  could  be  computed 
analytically  to  some  distance  R^.  Then  another  expansion 
could  be  made  about  R^,  and  a  whole  new  set  of  coefficients 
for  each  variable  simultaneously  computed.  The  new  expan¬ 
sion  could  be  checked  by  calculating  the  variable  back 
toward  the  axis  and  comparing  the  radial  profiles  of  the 
variables  expressed  by  both  series.  This  process  could  be 
extended  in  short  segments  to  any  radius.  A  possible 
limitation  would  be  a  region  where  all  the  variables  be¬ 
come  rapidly,  changing  functions  of  r,  in  which  case  the 
number  of  terms  in  the  Taylor  series  necessary  to  repre¬ 
sent  each  variable  accurately  might  be  prohibitive. 

Let  us  then  investigate  the  possibility  of  applying 
such  a  contrived  method  to  the  eight  differential  equa¬ 
tions.  Immediately,  the  equations  are  non-llneari  that  is, 
products  of  variables  appear  In  them.  So  the  simultaneous 
coefficient  equations  would  be  non-linear.  Although  the 
methods  of  solving  simultaneous  linear  equations  are  well 
established,  .the  generalized  methods  pertaining  to  non¬ 
linear  equations  become  exceedingly  complex  for  more  than 
Just  two  or  three  equations. 

However,  the  Identical  problem  arises  which  prevented 
the  ohange  of  variables  to  n  and  in  the  last  section. 
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Note  that  the  coefficients  of  the  n  and  ^  series  would 
still  have  to  come  from  the  simultaneous  solution  of 
Eqs  (56)  and  (57).  It  was  found  before  when  solving  for 
d^/dr  that  no  cancellation  of  terms  In  Eqs  (56)  and  (57) 
could  be  accomplished.  Hence  the  large  common  factor 
between  rig  and  n^  could  not  be  eliminated.  This  same 
problem  would  appear  In  the  simultaneous  solution  for  the 
series  coefficients  of  .  Terms  like  (E  +  *0  B)  and 
(E  +  UpB)  would  only  be  translated  Into  terms  Involving 
the  coefficients  of  E,  Vg  ,  Ug ,  and  B.  The  difficulty  of 
computing  relatively  small  differences  between  large  num¬ 
bers  would  remain.  This  dilemma  was  actually  verified  by 
expansion  of  all  variables  as  described  above. 

In  this,  section  It  has  been  shown  how  even  small 
errors  In  the  variables  make  the  numerical  integration  of 
the  equations  by  such  methods  as  Runge-Kutta  and  Milne 
Impossible.  The  Taylor  method  has  been  proposed  and  it  has 
been  demonstrated  not  feasible  also. 

Summary 

This  chapter  deals  essentially  with  the  failure  to 
find  a  method  of  accurately  integrating  the  equations 
derived  in  Chapter  III.  A  set  of  variables  Is  first  chosen 
from  which  the  variables  of  Interest  can  be  algebraically 
computed.  The  sources  of  error  are  identified  in  the  terms 
of  the  differentials  of  these  variables  so  that  the  reader 
can  understand  the  purpose  of  the  subsequent  arguments. 
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The  differential  equations  are  written  for  r=0,  and 
expressions  for  Initial  conditions  on  all  variables  are 
derived  from  them.  Such  solutions  are  required  before 
numerical  integration  methods  can  be  applied.  The  Initial 
condition  equations  are  used,  with  some  typical  values  of 
the  specified  constants,  to  illustrate  the  severity  of  a 
major  obstacle  to  the  solution  of  the  equations  away  from 
the  axis.  It  is  found  that  a  property  of  the  system  of 
equation  is  the  very  small  relative  difference  on  axis 
between  ng  and  Zn^ . 

The  results  of  the  application  of  several  common  inte¬ 
gration  methods  are  reported.  The  behavior  of  the  variables 
under  the  Milne  iterative  Integration  scheme  confirms  that 
some  of  the  differentials  are  subject  to  a  severe  loss  of 
accuracy.  The  electric  field  equation  is  most  immediately 
affected  by  the  loss  of  accuracy  in  the  difference  between 
the  particle  densities.  The  error  in  E  causes  the  differ¬ 
entials  of  the  density  terms  to  be  erroneously  computed. 

The  most  important  reason  for  this  is  the  equality  of  v  B 
with  (— E )  ,  and  u  B  with  (-E).  The  effect  of  error  in  the 
terms  is  amplified  by  the  cancellation  of  their  most 
significant  parts. 

No  rearrangement  of  the  differential  equations  could 
be  found  such  that  the  cancellation  of  the  common  factors 
between  nearly  equal  terms  can  be  performed  analytically. 

The  equations  came  originally  from  the  separate  application 
of  the  Boltzmann  equation  to  the  electron  and  ion 
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distributions.  It  is  this  separate  consideration  of  the  two 
types  of  particles  which  causes  the  dilemma  at  hand.  The 
two  distributions  cannot  be  handled  separately.  The  eleo- 
trlc  field  equation  and  the  force  term  Fc  unavoidably  tie 
the  two  together. 
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V.  Conclusion 

This  chapter  proposes  recommendations  for  continued 
study  In  the  areas  discussed  In  Chapters  II  through  IV. 
Plans  have  been  already  made  at  the  Plasma  Physics  Research 
Laboratory  of  the  Aerospace  Research  Laboratories  to  formu¬ 
late  methods  for  overcoming  the  difficulties  explained  In 
Chapter  IV,  and  these  plans  are  reported. 

Dlsplaoed-Coll  C onf  1 gurat 1 on 

T^e  study  of  the  magnetic  field  of  the  dlsplaced-coll 
configuration  as  reported  In  Chapter  II  could  be  greatly 
expanded.  It  Is  limited  In  validity  by  the  approximations 
of  Inf lnltely- thin  colls  In  an  infinitely  extended,  straight 
arrangement.  It  would  be  useful  to  devise  computer  programs 
similar  to  those  described  In  Appendices  B  and  C  to  map  the 
magnetic  field  of  a  system  of  real  solenoids  arranged  as  In 
Fig.  1.  The  general  solenoldal  magnetic  coll  can  be 
treated  mathematically  just  as  the  single  loop  of  current. 
Accurate  and  convenient  computer  programs  have  been  written 
to  oompute  the  magnetic  field  of  round  coils  anywhere  out¬ 
side  the  current  windings  (Refs  3  and  9). 

Besides  the  additional  parameters  Introduced  In  using 
real  oolls,  new  variables  would  arise  due  to  the  toroidal 
arrange  me  nt .  The  curvature  of  the  axis  about  which  the 
oolls  are  plaoed  would  obviously  destroy  the  azimuthal 
symmetry  of  the  magnetic  field.  The  gradients  In  field 


68 


GSP/PH/69-7 


strength  transverse  to  the  general  direction  of  the  mag¬ 
netic  field  should  increase.  The  most  profound  effeot 
would  probably  be  on  the  cyclic  rotation  of  the  B  lines. 

In  Chapter  II  it  was  demonstrated  that  such  charac¬ 
teristics  as  longitudinal  and  transverse  gradients  of 
magnetic  strength,  the  ratio  of  transverse  to  axial  B 
vector  components,  and  the  rotation  of  B  lines  are  func¬ 
tions  of  the  spatial  parameters  of  the  coil  configuration 
of  Pig.  2.  Very  broad  conclusions  were  drawn  about  the 
feasibility  of  selecting  specific  field  characteristics  by 
adjusting  these  spatial  parameters.  However,  should  a 
toroidal  device  such  as  that  depicted  in  Fig.  1  be  found 
promising  for  plasma  containment,  some  procedure  will  be 
needed  to  design  magnetic  fields  with  specified  charac¬ 
teristics  . 

The  choice  of  exact  sizes  and  shapes  of  solenoids  are 
somewhat  limited  by  what  is  commercially  available.  The 
designer  would  then  have  to  rely  primarily  on  the  spatial 
arrangement  of  the  colls.  Computer  techniques  probably 
would  be  the  best  tools  for  finding  directly  what  configu¬ 
ration  of  available  solenoids  would  yield  a  magnetic  field 
most  nearly  like  that  desired.  Many  of  the  computation 
techniques  developed  in  Appendices  A,  B,  and  C  would  be 
adaptable  in  such  procedures.  Another  thesis  done  in 
cooperation  with  the  Plasma  Physios  Research  Laboratory 
reports  the  attempts  of  Capt.  D.  B.  Taylor  to  design  ooil 
configurations  generating  speoiflo  magnetic  fields  by 
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using  standard  coils  (Ref  17il). 

Indirect  Solution  of  Differential  Equations 

It  became  evident  in  Chapter  IV  that  essentially  only 
one  difficulty  prevents  the  direct  numerical  integration  of 
the  differential  equations  derived  in  Chapter  III.  It  is  a 
characteristic  of  the  nearly  equilibrium  plasma  that  the 
relative  difference  between  ion  and  electron  number  densi¬ 
ties  is  extremely  small.  No  change  of  variables  could  be 
found  to  eliminate  the  necessity  of  subtracting  the  nearly 
equal  densities.  The  loss  of  accuracy  in  the  differences 
is  reflected  in  the  erroneous  calculation  of  important 
differentials . 

Lt.  Col.  Wingers on  has  proposed  an  alternative  to 
solving  the  equations  simultaneously  from  initial  conditions 
only.  It  Involves  a  mating  of  theoretical  prediction  and 
experimental  verification.  Instead  of  trying  to  compute 
the  difference  between  ion  and  electron  densities,  it  Is 
proposed  that  the  difference  be  Imposed  on  the  equations  as 
a  radial  function.  This  would  be  equivalent  to  specifying 
the  radial  electric  field.  On  the  basis  of  this  function 
and  the  arbitrary  injection  of  charged  particles  the 
system  of  differential  equations  might  be  solved  for  the 
variables  of  interest. 

These  variables  oan  also  be  indirectly  measured  in 
real  plasmas.  Plasma  probe  technology  is  suoh  that  parti- 
ole  densities  and  rotational  velocities  are  measurable  to 
at  least  an  order  of  magnitude  (Ref  6il).  From  a 
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knowledge  of  the  radial  density  profiles  the  continuity 
equations  can  be  used  to  compute  the  radial  drift  veloci¬ 
ties.  Although  the  magnetio  strength  In  a  hot  plasma  cannot 
be  measured  directly,  the  difference  between  the  field  due 
to  a  set  of  magnetic  colls  In  a  vacuum  and  In  the  presence 
of  a  plasma  can  be  predicted  (Ref  12»73).  So  the  strength 
of  the  magnetic  field  associated  with  a  plasma  can  be  cal¬ 
culated  Indirectly.  Temperatures  for  both  electrons  and 
Ions  are  also  available  by  probe  measurement. 

Thus  all  the  variables  which  appear  In  the  differential 
equations  can  be  measured  In  a  real  plaBma.  Careful  com¬ 
parison  of  the  mathematical  solutions  and  the  measurements 
could  show  whether  the  electric  field  Imposed  mathemati¬ 
cally  Is  consistent  with  reality.  The  most  Important 
Indication  would  be  whether  the  shape  of  the  radial  electric 
field  profile  Is  reasonable,  rather  than  whether  the  abso¬ 
lute  value  of  the  field  at  any  point  Is  exact.  The  contour 
of  the  electric  field  would  Indicate  what  type  of  series 
expansion  would  represent  the  field  most  accurately.  Then 
as  suggested  In  Chapter  IV  all  other  variables  could  be 
expanded  In  series  consistent  with  the  electric  field 
function,  and  the  entire  set  of  equations  solved  simultane¬ 
ously  for  the  coefficients  of  each  variable's  expansion. 

The  Plasma  Physics  Research  Laboratory  has  constructed 
an  axisymmetric  plasma  containment  device  called  the  ELMAX. 
ELMAX  is  equipped  with  Langmuir  probes  and  the  associated 
eleotronlcs  to  measure  radial  profiles  of  plasma 
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temperatures,  number  densities,  and  rotational  frequencies. 
Lt.  Col.  Wingerson  plans  to  use  the  ELMAX  to  obtain  the 
measured  variables.  The  configuration  of  the  device  is 
similar  to  the  system  hypothesized  In  this  paperi  it  has  an 
axial  magnetic  field,  and  gases  can  be  added  to  the  running 
plasma  along  the  axis  of  the  system.  Thus  the  equations  of 
Chapter  III  should  describe  macroscopic  phenomena  in  ELMAX 
reasonably  well. 
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Appendix  A 

Equations  of  the  Dlsplaced-Coll  C onf lgurat  1  on 


Thin  Loop  Eo  uat  1  ons 

The  radial  and  axial  components  of  the  B  field  of  an 
Inf lnltely- thin  loop  of  current  In  cylindrical  coordinates 
are 


Br  =  AoZ 


2.1 r  r  [[cl-*- r]  x  +  z]  1/aL  [cl- r] 


-k-t-lo? 


iN-rv+z*]  F 

_-rlx4-z.x  J 


( 9*+) 


]7  =  ^1  1 

k  +  03-  rx-z.a  F 

Xvr  [[a.+  r]x  +  z.xj 

[a.-r]  x+z.x  _ 

(95) 


(Ref  20il55).,  where  the  axis  of  the  loop  Is  the  Z  axis  of 
the  coordinate  system.  The  azimuthal  component  of  the  field 
Is  zero  by  symmetry.  "I"  Is  the  current!  "a"  Is  the  radius 
of  the  loop  in  meters  1  and  K  and  E  are  the  complete  elliptic 
Integrals  of  the  first  and  second  kind,  respectively. 

Normalization.  Let  us  normalize  units  such  that 
4j£l/27fa  Is  equal  to  unity.  Then  B  Is  a  dimensionless  number. 
All  coordinates  are  measured  as  fractions  of  the  radius  of 
the  loop. 

Elliptic  Integrals .  The  elliptic  Integrals  K  and  E 
are  defined  by  the  relations 


o 


75 


GS P/PH/6 9- 7 


E  =  J [l--ft1sin1e]  Vi  d 


9 


(97) 


where  k,  called  the  modulus  of  the  integral,  Is  given  by 

'/x 


-fc  =  _ ^Q-r 

_  [a+  r]  a  +  z?  J 


(98) 


The  complementary  modulus  ok  is  defined  such  that 

CK1  4  KX  -  \  (99) 

The  elliptic  integrals  can  be  expanded  into  polynomial 
approximations  of  various  numbers  of  terms.  Two  expansions 

O 

which  are  within  2  x  10“  of  the  true  values  of  K  and  E  for 
any  given  k  are  given  below  (Hef  li591). 


k[c  -h]  -  [a«  4-  oljCk1  4-  aaCK4  4  a3a<fc  4-  cu+ck8] 

4  [k  4  LhCK1  4  ^0^4  4  b^CK8]  InO^k] 

E[c-fc]  =  [l  4C,CKa  4  CjCK^  4  C3CKfe  4C+C.K8] 
4[d,CK1  4  dxCK^  4  d^vc*’  4  d^CK&]  InD^J 


(100) 


(101) 


These  approximations  are  useful  because  of  the  ability  of 
the  digital  computer  to  handle  such  calculations  with  great 

speed. 

Singularities .  As  mentioned  above,  k  must  be  less 
than  one.  The  cylindrical  coordinates  corresponding  to  k 
equal  to  one  are  r  equal  to  one  and  z  equal  to  zero,  as  can 
be  seen  Immediately  from  Eq  (  98  ) .  These  are  the 
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coordinates  of  the  conducting  loop  Itself.  This  Is 
expected  since  any  conductor  Is  a  singularity  In  the  mag¬ 
netic  field.  Otherwise  the  equations  for  Br  and  Bz  are 
valid  for  any  coordinates  r  and  z,  r  not  equal  to  zero. 


Superposition  of  Coll  Fie Ids 

The  principle  of  superposition  is  exercised  to  find 
the  magnetic  field  at  any  point  due  to  the  fields  of  a 
number  of  separate  colls.  Fields  of  several  colls  can  be 
considered  together  at  a  point  by  adding  vectorally  the  B 
vector  from  each.  This  Is  most  conveniently  accomplished 
by  adding  parallel  components.  An  orthogonal  coordinate 
system  Is  defined  In  Chapter  II  relative  to  the  coll  con¬ 
figuration. 

Position  Re lat Ive  t o  a  Coll.  The  axial  dlstanoe  from 
any  point  to  the  plane  of  each  coll  Is  easily  found.  It  Is 
simply  the  Z  coordinate  of  the  point  relative  to  the  system 
axes  plus  or  minus  an  Integral  multiple  of  the  dlstanoe 
between  the  planes  of  the  colls.  The  X  and  Y  coordinates 
relative  to  an  individual  coll  depend  on  the  position  of  the 
center  of  that  coll  with  respect  to  the  system  axes.  If 
the  transverse  coordinates  of  a  point  are  X  and  Y,  the 
coordinates  relative  to  a  coll  whose  center  lies  at  "dlsp" 
off  the  axis,  and  at  an  angle  "alpha”  from  the  X  axis  are 


Xc=  X-disp  Cosoc 
Vc  =  Y-  disp  sin  c*. 


(102) 


77 


y 


GS P/PH/6 9- 7 


Hence,  the  radial  distance  from  the  center  of  the  col3  to 
the  point  is 

Rc  ,  [Xe1  +  Yc2]  <™» 

Thus  the  R  and  Z  distances  can  be  found  for  calculation  of 
the  magnitudes  Br  and  Bz  of  the  field  of  any  coil  at  a  given 
point . 

Summat 1 on  of  Components .  The  Bz  components  from  each 
coil  are  added  since  they  are  parallel.  The  Br  components 
must  be  divided  into  Bx  and  By  components  which  are  parallel 
in  direction  to  the  system  axes. 

The  directions  of  the  Br  from  the  coils  are  radial 
outward  from  the  centers  of  the  coils  to  the  point  under 
consideration.  The  Br  are  divided  via  multiplication  by 
direction  cosines i 

,:Bx =  Br 


(105) 

In  summary,  approximate  expressions  for  the  axial  and 
transverse  components  of  B  for  a  single  thin  loop  have  been 
derived.  A  coordinate  system  is  used  as  the  basis  for 
dividing  these  components  into  Bx,  By,  and  Bz  components. 
Thus  the  fields  of  separate  coils  are  added  veot orally  at  a 
point . 


Xc 


Rc. 


(104) 
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Appendix  B 
The  Subroutine  Field 

This  appendix  describes  a  subroutine  called  "Field" 
which  is  designed  to  provide  the  magnitude  and  the  direc¬ 
tion  of  the  magnetic  strength  vector  B  at  any  point  in  the 
field  of  the  displaced-coil  oonf iguratlon.  The  origins  of 
the  requirement  for  the  subroutine  are  detailed  In  Chap¬ 
ter  II.  The  objectives  and  criteria  for  design  are  stated 
here.  The  method  of  computation  is  explained,  and  operat¬ 
ing  Instructions  are  provided.  Finally,  the  effectiveness 
of  the  subroutine  is  evaluated. 

Fig.  8  is  a  flow  diagram  of  Field,  and  Fig.  9  Is  a 
listing  of  the  computer  code.  Table  III  defines  the 
Fortran  variable  names  used  in  Field, 

Objectives  and  Criteria 

Bm  and  Direction  Cosines .  In  Chapter  II  it  Is 
determined  that  RTA ,  the  ratio  of  transverse  to  axial 
vector  components,  and  the  uniformity  of  Bm,  the  magnitude 
of  B,  should  be  studied.  Als.o,  the  Integration  of  field 
lines  requires  the  direction  cosines  of  B.  It  is  the  objec¬ 
tive  of  Field  to  supply  Bm  and  the  three  direction  cosines 
of  B  at  any  point  in  the  generalized  displaoed-ooil 
o  onf  iguratlon’, 

RTA  Is  the  square  root  of  the  sum  of  the  squares  of 
the  transverse  direction  cosines.  The  uniformity  of  Bm  can 
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Pig.  8.  Flow  Diagram  of  Field. 
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S I BFTC  FIELD.  DECK 

.  SUBROUTINE  F I ELD ( CO • DCO.BM. SPACE. DI SP .ALPHA I 
DIMENSION  CO ( 3 )  • DCO ( 3  ) 

AL PHAR“ ALPHA/ 57  • 3 
TOL- .001 
GAMMA«0. 

C  YCLE“SPACE 
ENTRY  FUN( N.S.CO.DCO J 
C  TEST  AXIAL  COORDINATE  OF  POINT 
IF  <  CO( 3 ) -CYCLE  >  2.3.3 
C  IF  ANOTHER  LOOP  HAS  BEEN  PASSED.... 

3  CYCLE  -CYCLE+SPACE 
gamma*gamma+alphar 

2  C  HI “GAMMA 

PSI  -GAMMA-f  ALPHAR 

C  FIND  Z  COORDINATES  RELATIVE  TO  NEAREST  LOOPS 
ZLEF«CO<  31-CYCLE 
ZRIT-ZLEF-SPACE 
BX»0. 

B  Y“  0  • 

BZ«0. 

C  FIND  X-Y  COORDINATES  RELATIVE  TO  TWO  LOOPS 
IO  X  “CO< I )-DISP*COS(CHI > 

Y  -CO< 2I-DISP*SIN<CHI I 
R2“X*X-.Y*Y 
R-SQRT(R2) 

Z2-ZLEF*ZLEF 

T  OP  1“1«+R2+Z2 
TOP2- l 1.-R2-Z2) 

BOT I  * SORT  < <1.  +  R)*<1.+R>+Z2> 

BOT2-<  l.-R)*U.-R>  +  Z2 
AK»4.*R/C(1.+R)*C1.+R)+Z2I 
CALL  CE I  PA (AK.FEI ,SEI  ) 

BR1“(ZLEF/RL/B0T1>*  <  -F  E  I  +  TOP  1  /BOT2*SE  I  ) 
BZ1“<FEI+T  OP 2 /BOT  2*S  E I I/BOTl 
X  «CO< 1I-DISP*C0S<PSI  ) 

Y  “CO ( 2 > -D I SP*5 1 N { PS  I  ) 

R2“X*X«-Y*Y 

R-SQRTIR2) 

Z2-ZRI T*ZRIT 
T  OP 1“1.+R2  +  Z2 
T  OP  2 “ ( I.-R2-Z2) 

BOT  1  “SORT  (  <1.*-R)*(1.  +  R)+Z2> 

BOT2  “  (  l.-R)*»l.-R>+Z2 
AK»4.*R/< ( l.+R )* ( 1 .-*-R  }*Z2  ) 

CALL  CEIPACAK.FEI .SEI  > 

BR2“<ZRI  T/RR/BOT  I  J  *  (-FE1  *•  TOP  1 /B0T2*SE  I  ) 

BZ2“ • FE I +  TOP2 /BOT  2*S  E I )/BOTl 


Fig.  9.  Listing  of  Field 
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BX«BX-t-Xl/RL*BRl  +  X2/RR*BR2 
BY«BY+Y1/RL*BR1+Y2/RR*BR2 
BZ-BZ+BZI+BZ2 

C  TEST  SIZES  OF  AXIAL  CONTRIBUTIONS 
IF<BZl-TOL*BZ )  21.21,29 
IF(BZ2-TOL*8Z )  22,22,2 9 
C  IF  MORE  COILS  NEED  TO  BE  CONSIDERED.... 
29  ZLEF-ZLEF+SPACE 
ZR I T*ZR I T-SPACE 
CHI »CHI— ALPHA R 
PS I  * PS  I +ALPHAR 
GO  TO  10 
22  CONTINUE 

BM«SQRT<  BX*BX+BY*BY  +  BZ*BZ  » 

DCO< 1 ) ■BX/8M 
DCO ( 2  I »BY/BM 
DCO( 3 J-BZ/BM 
RETURN 
END 


Pis*  9*  Listing  of  Pield,  (oont.) 
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Table  III 


Definitions  of  Fortran  Variables  In  Field 
Fortran  Name  Definitions 


AK 

ALPHAR 
BX,  BY,  BZ 
BR1,  BR2 

BZ1,  BZ2 
CHI,  PSI 

CK 

CO 

CYCLE 

DCO 

FEI ,  SEI 

GAMMA 

N 

R2 ,  Z2 
RL,  RR 

S 

TOL 

TOPI,  TOP2, 
B0T1,  BOT2 


Modulus  of  elliptic  Integrals 

Alpha  In  radians 

Vector  components  of  B 

Transverse  terms  of  B  from  a  coll 
pair 

Axial  terms  of  B  from  a  coll  pair 

Azimuthal  coordinates  relative  to  a 
coll  pair 

Complementary  modulus  of  elliptic 
Integrals 

3-mat rlx  of  position  coordinates 

Integral  multiple  of  space 

3-matrlx  of  direction  cosines 

Elliptic  integrals,  first  and 
second  kind 

Integral  multiple  of  ALPHAR 

Number  of  dimensions  (3) 

Squared  radial  and  axial  distances 

Radial  distances  from  point  to  coil 
centers 

Integral  length  of  B  line 

Parameter  of  accuracy  required  of 
Bz  term 

Terms  In  the  equations  for  Br,  Bz 
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Definitions 
Fortran  Name 

X,  Y,  ZLEF ,  ZHIT 


Table  III  (cont.  ) 


of  Fortran  /ar tables  In  Field 
Definitions 


Field  point  coordinates  relative  to 
a  coll  pair 
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be  studied  by  computing  Bm  as  a  function  of  changing  radial 
and  axial  position.  For  convenience  and  generality,  there¬ 
fore,  Field  is  written  as  a  subroutine  to  be  called  by  any 
program  which  requires  the  magnitude  and  direction  of  B. 

Field  is  designed  to  provide  for  the  control  of  errors 
of  approximation.  The  necessity  for  error  control  arises 
primarily  from  the  possibility  that  in  an  integration 
process  using  the  output  of  Field  errors  may  be  cumulative. 
For  this  reason  the  maximum  errors  in  the  direction  cosines 
are  firmly  under  control. 

Method  of  Calculation 

This  section  explains  the  application  of  previously 
derived  equations  to  the  objectives  outlined  above.  It 
also  describes  the  mechanism  of  error  control. 

Vector  Component  Computation.  Equations  are  developed 
in  Appendix  A  for  the  normalized  components  of  B  in  the 
field  of  a  single  loop  of  current.  The  coordinates  of  both 
the  field  point  and  the  coll  center  must  be  known  in  order 
to  apply  these  equations. 

Adding  Components .  Field  Initially  determines  from 
the  coordinates  of  the  specified  point  the  coordinates  of 
the  centers  of  the  two  colls  which  bracket  the  point. 
Calculations  of  vector  components  are  oarried  out  on  one 
pair  of  coils  at  a  time,  beginning  with  the  pair  closest  to 
the  point.  Pairs  successively  farther  out  from  the  point 
in  both  axial  directions  are  considered  in  turn.  The  pairs 
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of  vector  components  are  added  to  the  subtotals  of  Bx,  By, 
and  Bz  as  they  are  calculated. 

Testing  Axial  Components .  The  Individual  axial  com¬ 
ponents  are  expected  to  decrease  with  Increasing  distance 
between  the  point  and  the  coll.  Field  tests  each  axial  - 
contribution  against  a  value  called  "tol.  ••  Tol  is  an  error 
parameter  defined  as  the  minimum  relative  size  of  axial 
contribution  to  be  considered.  It  Is  approximately  equal 
to  the  difference  between  the  total  axial  strength  at  a 
point  in  an  Infinite  system  of  coils  and  the  total  of 
oonti ibutlons  from  a  finite  number  of  coils. 

If  both  axial  terms  from  a  pair  of  coils  are  not 
smaller  than  tol.  Field  considers  the  outer  next  pair  of 
colls.  If  the  test  is  met,  however,  the  additive  process 
Is  terminated.  The  absolute  vector  magnitude  and  the 
direction  cosines  are  computed.  Control  is  returned  to  the 
calling  program. 

Operating  Instructions 

The  following  control  statement  must  be  placed  -in  any 
program  calling  Field  i 

DIMENSION  CO(3),DCO(3) 

where  CO  is  the  array  of  point  coordinates,  and  DCO  is  the 
array  of  direction  cosines.  Space,  disp,  and  alpha  must 
be  defined,  and  N  should  be  set  equal  to  the  Integer  three 
in  the  calling  program. 
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The  first  time  Field  is  called  for  a  new  combination 
of  geometric  parameters  the  following  statement  is  appro¬ 
priate  i 

CALL  FIELD(CO,DCO,BM, SPACE, DISP,  ALPHA) 

Bm  and  the  direction  cosines  are  available  in  the  argument 
list  upon  return  from  Field.  This  statement  must  be  used 
if  Bm  is  required  at  eaoh  point. 

If  Field  is  serving  a  line  integration  routine,  the 
following  statement  should  be  used  for  each  integration 
step  after  the  above  statement  has  been  used  once  i 
CALL  FUN(N,S,CO.DCO) 

Error  Control.  Tol  is  set  within  Field  to  a  value 
which  insures  a  maximum  error  in  Bm  of  one  tenth  of  one 
percent.  If  a  lesser  error  is  desired,  tol  can  be  reset 
within  Field i 

Evaluat 1 on  of  Effectiveness 

This  section  is  an  evaluation  of  how  well  Field  can  be 
expected  to  meet  the  objectives  and  criteria  previously  set 
down.  Any  stated  observations  are  based  on  experience  in 
using  Field. 

Limitations  on  accuracy  of  results  from  Field  will  be 
considered  first.  The  absolute  limitation  upon  accuracy 
arises  from  the  approximate  method  of  computing  the 
complete  elliptic  integrals.  It  is  noted  in  Appendix  A 
that  the  method  ohosen  has  a  maximum  error  of  2  x  10"®  in 
each  integral  calculation.  It  is  therefore  concluded  from 
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a  study  of  the  component  equations  of  Appendix  A  that  the 
minimum  meaningful  tol  that  should  be  considered  is  10-^. 

The  specification  of  tol  should  be  a  compromise 
between  desires  for  accuracy  and  economy  of  computer  time. 
It  is  advised  that  a  reduction  of  tol  by  one  order  of 
magnitude  roughly  doubles  the  number  of  coils  that  con¬ 
tribute  at  the  field  point  relative  axial  components 
greater  than  tol.  Note  that  with  the  method  used  by  Field 
the  calculations  of  the  contributions  of  a  single  coil 
require  the  same  number  of  operations  wherever  the  coil  is 
located  relative  to  the  specified  point. 

It  was  mentioned  earlier  that  tol  has  been  preset 
within  Field  such  that  only  coils  which  contribute  at  least 
0.1  percent  of  the  total  field  Strength  at  a  point  are  con¬ 
sidered.  This  does  not  mean  that  the  three  direction 
oosines  can  be  in  error  by  this  amount.  It  was  found  that 
the  consideration  of  a  coll  far  enough  away  axially  from  a 
point  to  contribute  less  than  C.l  percent  of  the  field  due 
to  all  coils  changed  any  single  direction  cosine  by  less 
than  10“^  peroent.  The  reason  for  this  is  that  the  fields 
due  to  colls  far  from  a  point  are  nearly  perfectly  axial  in 
direction.  The  transverse  component  from  a  single  coil  far 
away  is  very  small  so  that  it  adds  very  little  to  Bx  and 
By.  Also  sinoe  the  field  from  a  distant  coil  is  nearly 
axial  in  direction,  its  total  contributed  magnitude  is 
almost  entirely  due  to  the  axial  component.  Note  that  the 
axial  direction  cosine  is  Bz/Bm,  and  that  the  field  of  the 
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entire  dlsplaced-coil  configuration  at  any  radius  is  very 
nearly  social.  Then  Bz/Bm  is  nearly  equal  to  one,  and  Guid¬ 
ing  nearly  equal  amounts  to  Bz  and  Bn  changes  the  ratio 
very  little . 

Summary 

Field  is  capaoie  of  providing  to  a  calling  progreun  the 
magnitude  and  direction  cosines  of  B  at  any  non-singular 
point  in  the  magnetic  field  of  the  dlsplaced-ooll  configu¬ 
ration.  Error  control  is  accomplished  by  bringing  the 
axial  component  of  B  to  within  a  certain  percentage  of  the 
axial  strength  expected  in  an  infinite  system  of  oolls. 
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Appendix  C 
The  Pr ogram  Tracer 

This  appendix  describes  a  program  called  "Tracer" 
which  Is  designed  to  compute  the  paths  through  space  of 
field  lines  in  the  magnetic  field  of  the  dlsplaced-coll 
configuration.  The  appendix  reiterates  the  objectives  of 
the  study  for  which  Tracer  is  designed,  and  it  describes 
the  methods  by  which  Tracer  meets  those  objectives. 

Operating  instructions  ere  provided  for  using  Tracer.  The 
appendix  is  concluded  with  an  analysis  of  the  control  of 
errors  of  integration.  Fig.  10  Is  a  flow  diagram  of  Tracer, 
and  Fig.  11  is  a  listing  of  the  computer  code. 

Objectives  of  Tracer 

It  is  very  useful  to  Investigate  two  properties  of  the 
magnetic  field  of  the  displaced-coil  configuration.  First 
the  field  lines  are  expected  to  exhibit  a  pronounced  helical 
twist  similar  to  the  helix  on  which  the  coil  centers  lie. 
This  characteristic  is  called  the  first  order  twist.  Also 
secondary  rotations  of  the  entire  pattern  of  lines  about 
the  axis  may  exist. 

There  then  is  a  need  for  a  computer  program  to  numer¬ 
ically  integrate  the  paths  of  the  magnetic  field  lines  of 
the  coll  system.  It  is  the  objective  of  Tracer  to  inte¬ 
grate  the  paths  of  field  lines  which  begin  at  selected 
coordinates  in  the  zero  plane  of  a  parametrically 
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Fig.  10.  Flow  Diagram  of  Tracer 
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* 


$  I  BF  TC  TRACER 

DIMENSION  XSTRTl  1C!  .  YS  T:<T  (  1  ?  J  .CO  <  3  > 

EXTERNAL  FUN 
LOGICAL  VCDL 
REAL  LIVIT 

1  READ!  5.81)  SPACE. DI  ,L  .C* . 

1  (XSTRTH  )  .wrTRTI  I  J  •  I'!  .1''  J 
81  FORMA  T(6X»2FF.2»Ft.1.  .  !'».F?  .<?/ 
15X.10F?.2/fX,l''Fc.;; 

WRITE  (6*92)  SPACE.-’  • 1  -  «  L 

92  FORMA TC 1H1.1CX.19  :?ACC  .  T  .  r  ,  * 

110X.  2  7MDI c  Pi  A  CEF  •  '  r  ?«  COIL  CL  :T;  . 

2 1 ox . 2 1 h anc l c  '-:t„'.  \  :  .*u::  ,  ..  , 

310X.  ISHN'JMSER  OF  L  I  *!'  C  *  f  A /// * 

ALPHAR=ALP  A/37.3 

LIvIT*CN*Er,ACE*3bO./*LPni 

DI S=.OO0C1 

DSVIN*.05 

DSMAX=1.C 

N  =  3 

DO  13*0  I  1  =  1. L 
WRITE (6 *99) 

19X.6HLENGTM.  AX.AHSI7E.  10X  ♦  l'  X  .•  x  .  1  <v  7 

s=o. 

DS-.20 


CO  C 1 ) *XSTRT<  I  I  ) 

CO( 2 ) = YSTRT  t  : I ) 

COC3J*C. 

call  FICLDCCOtPCOt  "•5PAC  .  :  • • 

C  CALL  SET....I*;IT!'TES  :  ‘TEC"  '• T  l  ?V 

call  set r:,: ,:c.o.  n.o: - 

C  A  SAMPLE  OUTPUT  FORMAT.. ■.LE-TT'*.  *■  TCP 

20  MRITC(4«91) S»DS«COf 1 ) »COI 2) »LC(  1 

91  FORMAT (5X.EF1C.A) 

c  call  step....I':tecr^t:s  _  _  :.r 

CALL  STEP 

21  IFCCOm-LIMIT) 

WRITE! 6»9C) 

90  FORMA  Til  HI) 

GO  TO  20 


29  -CONTrNUE 
100  CONTINUE 


4*  ) 

. 

'  .r, 


r 


Pig.  11.  Listing  of  Tracer 
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Table  IV 


Definitions  of  Fortran  Variables  In  Tracer 
Fortran  Name  Definitions 


CN 

DCO 

DIS 

DSMAX  (DSMIN) 
DS 
J 
L 

LIMIT 

MODE 


N 

S 

XSTRT  ,  YSTRT 


Number  of  first-order  cycles  covered 

3-matrix  of  direction  cosines 

Parameter  specifying  accuracy  In  STEP 

Maximum  (minimum)  size  of  DS 

Increment  of  Independent  variable 

Number  of  Integrating  steps 

Number  of  field  lines  to  be  Integrated 

Maximum  axial  distance  of  Integration 

Logloal  constant  specifying  variable 
(fixed)  DS  If  .False.  (.True.) 

Number  of  differential  equations  (three) 

Size  of  Independent  variable 

X-Y  coordinates  of  starting  point  of 
field  line 
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adjustable  eoll  configuration.  It  is  intended  that  the 
user  be  able  to  control  both  the  lengths  of  the  lines  and 
the  maximum  computational  error  that  can  arise  in  the 
coordinates  of  points  on  the  lines. 

Method 

A  magnetic  field  line  is  by  definition  in  the  same 
direction  as  the  local  magnetic  strength  vector.  Field 
lines  may  be  generated  mathematically  by  successively  inte¬ 
grating  over  an  Incremental  distance  the  direction  cosines 
of  the  local  B  vector. 

The  Subroutine  DFEQ.  Tracer  calls  upon  a  subroutine 
called  "DFEQ"  written  by  Mr.  Paul  J.  Nikolai  of  the  Aero¬ 
space  Research  Laboratories  at  Wright- Patters on  Air  Force 
Base,  Ohio.  The  subroutine  is  written  for  the  Fortran  IV 
computer  language. 

DFEQ  is  designed  to  Integrate  in  a  step-wise  manner  a 
set  of  simultaneous  first-order  differential  equations.  It 
requires  a  subroutine  written  by  the  user  to  supply  the 
differentials  of  the  dependent  variables.  Field  is  the 
necessary  subroutine  for  this  study.  Consult  Appendix  B 
for  a  description  of  Field. 

DFEQ  Integrates  the  first  three  steps  of  a  series  using 
a  classical  Runge-Kutta  method  (Ref  15«56).  A  four-point 
Adams-Bashf  orth-Adams-Moulton  predict  or- correct  or  scheme  is 
applied  to  the  fourth  and  succeeding  points. 

DFEQ  oan  be  used  in  a  very  useful  mode  to  vary  the 
size  of  the  integrating  increment  in  order  to  meet  a 
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certain  tolerance  for  accuracy.  Denote  by  YP  and  Y°  the 
predicted  and  corrected  values,  respectively,  of  the  1th 
dependent  variable  at  any  given  point.  "Dls"  1b  a  parameter 
specified  In  the  program  calling  DFEQ.  Define 

[y.c-  y»' 


err 


jua  , 

1 4-  moLX  [  l 


i,  7. 


(106) 


If  err  Is  greater  than  dls  for  any  *1*  the  step  size 
currently  used  by  DFEQ  is  halved,  and  the  test  is  applied 
again  over  the  same  region.  If  100  x  err  Is  less  than  dls, 
the  step  size  Is  doubled  for  succeeding  integrations.  Hence, 
DFEQ  Incorporates  a  mechanism  for  control  of  the  error 
Inherent  in  the  approximate  Integration  method  It  employs. 

Dls  Is  preset  for  the  convenience  of  the  user  to  a 
value  which  guarantees  an  agreement  of  .001  percent  between 
the  predicted  and . corrected  values  of  any  position  coordi¬ 
nate.  This  tolerance  Is  more  than  adequate  to  prevent 
cumulative  error  great  enough  to  obscure  the  actual  paths 
of  the  lines. 


Operating  Instructions 

This  section  gives  Instructions  necessary  to  the  use 
of  Tracer.  Field  must  be  complied  in  the  computer  along 
with  Tracer.  The  read  statement  in  Fig.  11  relates  how  the 
data  cards  must  be  prepared.  "L"  is  the  number  of  lines  to 
be  lntegratedi  up  to  ten  lines,  starting  at  the  coordinates 
(XSTRT ,YSTRT )  in  a  zero  plane,  may  be  traced.  The  distance 
over  which  the  lines  are  to  be  traoed  is  specified  by 
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selecting  ft  multiple  number  'CN)  of  cycles  of  coll  dis¬ 
placement  . 

Fig.  6  Is  an  example  cf  the  output  listing  provided  by 
Tracer.  It  is  a  simple  tabulation  of  the  progress  of  the 
Integration  for  a  single  line. 

Comments  on  Accuracy 

The  section  on  method  explains  that  dls  Is  specified 
within  Tracer  to  provide  for  a  maximum  relative  error  per 
step  of  10“ 5  in  the  coordinates  of  points  on  the  field 
lines.  On  the  basis  of  statistical  random  error  after  say 
N  steps  the  relative  error  in  any  field  line  coordinate 
should  be  less  than 

tcT5  Vn~ 

From  studies  of  the  average  drifts  of  field  lines  over  one 
hundred  cycles  of  first  order  twist  fcr  various  coil  con¬ 
figurations  it  was  determined  that  the  error  given  above 
amounts  to  less  than  . CC1  of  the  average  magnitude  cf  the 
drifts.  Since  only  a  qualitative  understanding  of  the  slow 
drifts  was  sought,  dls  was  not  reduced  to  attain  greater 
accuracy  in  the  integration  of  field  lines.  For  the  study 
of  first  order  field  line  twist  a  relative  accuracy  of  10“5 
per  step  was  Judged  more  than  satisfactory. 


Summary 


In  summary,  Tracer  is  a  program  which  integrates  the 
paths  of  field  lines  in  the  dlsplaced-coil  configuration. 
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The  method  of  calculation  provides  for  the  control  of  error. 
The  configuration  of  colls,  the  number  of  lines,  their 
lengths,  and  starting  points  are  selected  by  the  user. 
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Appendix  D 


Collision  Force  Between  Distributions 


Mention  Is  made  in  Chapter  III  of  a  term  which  In  the 
MHD  equations  can  be  used  to  account  for  Interactions 
between  particle  distributions.  The  nature  of  this  inter¬ 
action  Is  that  of  a  drag  force  arising  from  collisions 
between  different  particles.  Treated  as  an  average  the 
force  between  distributions  can  be  found  as  a  moment  of 
either  distribution  function.  Define 


(Ref  16i157),  where  df/dt  Is  the  change  due  to  collisions, 
and  the  integration  must  be  over  all  possible  values  of 
particle  velocity.  This  appendix  deals  with  the  evaluation 
of  this  Integral. 


Collision  Cross  Section 

Define  (v  v' )  as  the  geometric  cross  section  for  a 
change  in  particle  velocity  from  v  to  v*  due  to  an  encounter 
with  a  different  kind  of  particle.  Assuming  a  simple 
Coulomb  collision  between  the  particle  and  another  of 
veloolty  u,  the  cross  section  In  center-of-mass  coordinates 
Is 


q  sin  ^ 


(108) 
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(Ref  I81I50).  ©  is  the  scatter  angle  measured,  between  v 

and  v* ,  and 

=  [ul  — v]  ao9) 

q  and  q,  are  the  charges  of  the  struck  and  striking  parti¬ 
cles,  respectively,  and  M  is  the  reduced  mass  of  the  system 
of  the  two  particles. 


Force  Per  Particle 

Call  f(v)  and  F(u)  the  local  distribution  functions  of 
the  struck  and  striking  particles,  respectively.  The  flux 
of  bombarding  particles  relative  to  a  system  of  particles 
traveling  at  v  is  then  £*F(u).  Using  these  definitions  the 
probability  per  unit  time  of  v  changing  to  v*  due  to  a  col¬ 
lision  is 

P[v-»v']  =  g  F(uJ  e(c^,e)  (110) 

Define  the  change  in  momentum  per  unit  time. 

d  =  m  [v'-  v]  P  (v  — »  v'j  (m) 

It  is  necessary  in  the  integration  of  Eq  (111)  to 
account  for  all  cases,  represented  by  the  possible  values 
of  v*.  Assume  for  a  moment  that  v  and  u  are  known  before 
a  collision.  The  magnitude  of  v'  is  restricted  by  the  con¬ 
servation  of  energy  and  momentum!  it  is  in  principle  known 
from  v  and  u.  Consideration  of  all  possible  directions  of 
v* ,  or  equivalently,  all  possible  ©  would  therefore  be 
tantamount  to  consideration  of  all  v* .  We  therefore 
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Integrate  over  the  possible  scatter  angles  0  to  evaluate 

<%>• 

It  is  convenient  to  transform  to  a  center- of -mass 
coordinate  system.  Define  the  following t 
m  ■  mass  of  the  target  particle 
m0  ■  mass  of  impending  particle 
V  .*  (mv  +  me u)/(m  +  m0) 

M  =*  m  m0/(m  +  mQ ) 

Then  becomes  , 


<^>=  F(iJ  d 


©  (112) 


Define  also  an  orthogonal,  spherical  polar  coordinate 
system  with  £  as  one  of  the  orthogonal  axes  such  that 

=  3[^ CosQ  +^sin0Co5  cp  +  n  sine  sirv  £pj  -  (113^ 

This  definition  merely  puts  (£*  -  jg)  in  a  more  useful  form. 

(3'3]*a3siAf  ['Ssm-S 


(114) 


■+*  5  C.O.S  -S-  Cos  if  ft  Cos-S.  sin  <pj 

Now  define^  F(£,V)^  as  the  average  force  per  particle 
due  to  oollislons.  It  can  be  evaluated  as  the  volume 
integral  of 


;  <F(2-V)>=  /a<i>dn 


(115) 


100 


GS  p/m/6  9-  7 


where 

d_fl  =  2.  sin  Cos  -g-  d©  d<p  (116) 


is  the  differential  solid  angle  in  a  spherical  coordinate 
system.  Only  the  §  Integral  is  non-zero. 


<tkvj> 


A 


VnMg  [_  £o 


a  F(a)  in  [sin-S.  min]  (n7> 


where  is  the  scatter  angle  corresponding  to  the  so- 

called  cut-off  distance  for  collisions.  Define 


yv 


min 


(118) 


Spitzer  discusses  the  arguments  concerning  the  choice  of  the 
cut-off  distance.  His  resulting  equation  for->\.in  mks 
units  is 


_A_  »  12-n-  r^3k3T3lVl 

<£%  L  "  -I 


cnso 


(Ref  16 i 127). 


Volumetric  Force 

Eq  (117)  then  gives  the  time-averaged  force  per  parti¬ 
cle  due  to  Coulomb  collisions  between  unlike  particles. 
Define  Fc  as  the  volumetric  average  force.  It  can  be 
evaluated  as  the  moment  of  the  target  particle  distribution. 


Fc  =  |-p(y]  F|q,V)  a3Vd3g 


(120) 
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In  order  to  calculate  the  components  of  this  force  In 
cylindrical  coordinates  the  integral  Is  dotted  with  the 
unit  vector  *2“  In  the  direction  of  the  desired  force  compo¬ 
nent  . 


(121) 


Before  attempting  any  integrations  It  is  use  1  to 
choose  V  in  cartesian  coordinates  and  £  in  spherical  coor¬ 
dinates.  Then  the  respective  differential  volumes  are 


dV  =  dVr  dVe  dVz 


dcj  =  sm  ip  d©  d  cp  dg 


(123) 


The  limits  on  the  integration  of  these  differentials  are 

-oo  /  Vr,  V©/  Vz  +  oo 

v 

o<  it 


ox<g<  °° 

In  cylindrical  coordinates  the  three  dot  products 

•  r  =  sin  if  Cos  © 

•  ©  =  sin  cp  sin  © 
g  •  £  =  Cos  if 


3 


: h  g  are 

(124) 

(125) 

(126) 


102 


GSP/PH/69-7 


The  generalized  form  of  the  distribution  functions 
which  apply  under  the  assumptions  in  force  has  been  derived 
in  Chapter  III.  The  target  particle  distribution  is  given 
as  an  example . 


n 


'By 


m 


Lair  J 


exp 


-A.  /&rc\[\L 
2 


+  r  xii»]a 


(127) 


Recall  the  definitions  of  the  center- of -mass  velocity  V  and 
the  relative  velocity  The  product  f(v)F(u)  in  terms  of 

these  is 


V-  [/?-/£]  Mg  4-^cum  r  &  j  ^  {128) 


-  JL  |"q  -  fu)0-u)]  r  ©]al 

i  [Bn+A™ J  u  L  J  J] 

Integration  of  Fc^,  over  the  three  components  of  V  yields  a 
factor  of 


2.TX 


13/x 


l/Sm  +A  mQ  J 


(Ref  8,284). 

To  simplify  the  remaining  Integral  in  £  define  the 
following. 


J- 


3 


/S/Sc>rr\Yr\a _ 

_2  /£rr\  +  /5om0_ 


(129) 
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/o  =  r  [i 


j] 

/?rr\  + 


(130) 


Then  the  second  exponential  term  in  the  product  f(v)F(u) 


can  be  rewritten. 


-if  q  -  r  fuOo-Uj]  ©1 

2.  j_/5Vr\  +  _  J 

=  -  ~  2.  ^sin  ©  S  in  <p  -f  /?  aJ 

=  -  j[^-  /£ sin©  sin  <p]a4-  /?a  [1-  sina©  Sin\p 


(131) 


(132) 


The  integral  for  Fcm  is  then 


>A  nr>0  f~m  +mc 


itm  =•  f ^H-O  2  In-A.  nr>0  fm  +mo] 

L^-TTfcoJ  [/$m  +  /^omo] 

'  Jd©  Jd<p[g*m]  sin  tp  S»n©  sin  tp]  ^  (133) 

4-  2  [l  -  Sin2-©  Sin^l 

Consider  separately  the  three  integrals  corresponding 
to  the  three  different  terms  of  ^  •  'm.  The  'r  and  "z"  dot 
products  integrate  to  zero  by  properties  of  odd  functions 


of  &  and  .  The  ©  component  Integral  is 


X©  =  j~ d©  y"d(p  Sin^cp  Sin©  (134) 

rfi,sinesin<p  r  n 

-J  dJ7  exp- [  2  +  sina©  si  nacp]  J  d35) 
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el  Pe  ^  dx  -  JL 

a  I 


=  2rr 

Po 

\s  ^ 

Iq  can  be  written  in  terms  of  the  error  function,  ^ 

Ie=  -  /?  $'  (/?!] 

(Ref  10 i 696).  Although  £  is  not  analytic,  series  expan¬ 
sions  can  be  used  to  approximate  it  if  the  argument  ^  is 
less  than  one. 

Z__.  nl  [Hn  +1 


(136) 


(137) 


X  IT 1/1 


(138) 


n-o 

Therefore  the  complete  expression  for  the  only  finite  com¬ 
ponent  of  Fc  is 


Fc  =  © 

rav 

a  ln_A.nno 

mo 

-  - 

Xtt  ^ 

[^rn  -4-  Porr\o] 

30 


n  +  A  ^pXn-iL 


(139) 


Li) _ _ 

jpl.n-b.lJ  [n-1  l 


In  summary  a  volumetric  average  force  has  been  derived 
which  can  be  used  to  account  for  the  inter-distribution 
force  arising  from  Coulomb  collisions.  This  force  turns 
out  to  be  in  the  azimuthal  direction  only,  and  from  the 
definition  it  depends  on  the  difference  between  the 
azimuthal  drift  velocities  of  the  two  particle  distributions. 
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